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PREFACE 

This  Note  records  the  results  of  a  review  and  analysis  of  the  status  of  the 
computational  fluid  dynamic  (CFD)  modeling  techniques  related  to  the 
National  Aerospace  Plane  (NASP)  operation.  It  was  rmdertaken  as  a  task  in 
the  study,  "The  National  Aerospace  Plane  (NASP);  Development  Issues  for 
Follow-on  Systems,"  performed  within  the  Technology  Applications 
Program  of  Project  AIR  FORCE.  The  research  sought  to  evaluate 
independently  the  degree  of  uncertainty  and  the  technical  risk  involved  in 
predicting  the  NASP  performance  using  numerical  simulation  of  aerothermal 
and  chemical/combustion  processes.  This  Note  covers  the  technical  review 
portion  and  identifies  the  areas  for  research  emphasis  so  that  the  predictive 
reliability  of  the  NASP’s  potential  performance  parameters  can  be  improved. 

Sponsored  by  the  Air  Force  Directorate  of  Program  Planning  and 
Integration  (SAF/AQX),  the  project  has  also  required  close  cooperation  with 
the  NASP  Joint  Project  Office  (JPO). 

This  study  should  interest  those  concerned  with  the  aerospace  plane 
development  in  general,  and  hypersonic  CFD  modeling  in  particular. 


SUMMARY 


Computational  fluid  dynamic  (CFD)  simulations  are  now  gradually 
replacing  many  ground  experiments,  becoming  one  of  the  major  design  tools 
for  aerospace  engineers.  There  are  several  reasons  for  using  CFD  models; 

•  Technical:  Computational  limits  on  speed  and  memory  are 

rapidly  decreasing  with  time,  but  the  limits  on  experimental 
facibties  (wall  effects,  distortion,  etc.)  have  not  decreased. 
Economical:  Computer  speed  has  increased  faster  than  its  cost. 

Energy  saving:  Wind  tunnel  experiments  consume  vast  amounts  of 

energy, 

•  Convenience:  Computational  results  are  immediately  obtainable 

but  experimental  results  are  difficult  to  measure,  calibrate,  and 
interpret. 

Potential  application  and  major  hypersonic  aerodynamic  issues  that  may 
be  addressed  by  CFD  modeling  are: 

•  L/D  (lift/drag)  ratio  of  a  given  airframe  configuration, 

•  The  aerodynamic  heat  distribution  vs  maneuverability, 

•  Performance  parameters  of  the  supersonic  ramjet  (SCRAMJET), 

•  Aerodynamic  stability  within  the  NASP  flight  envelope, 

•  Boundary  layer  transition, 

•  Boundary  layer/shock  interaction, 

•  Other  viscous  effects. 

For  transonic  or  low  Mach  number  aircraft  design,  engine  and  airframe 
simulations  are  usually  made  separately.  But  for  the  hypersonic  plane,  the 
entire  aircraft  has  to  be  schematized  into  a  single  model,  because  the  highly 
integrated  hypersonic  plane  design  requires  dynamic  coupling  of  the  engine 
and  the  airframe.  A  model  containing  the  entire  aircraft  and  the  internal 
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combustion  computation  process  in  a  SCRAMJET  engine  would  eliminate 
problems  associated  with  the  specification  of  redundant  model  boundary 
conditions.  Consequently,  there  is  a  tradeoff  between  a  lower  spatial 
resolution  for  both  airframe  and  the  SCRAMJET  simulation  and  a  higher 
level  of  predictive  uncertainty  for  a  given  computational  resource. 

To  construct  a  CFD  model  covering  a  complete  hypersonic  plane  with 
reasonable  grid-space  resolution,  the  model  needs  substantial  computing 
resources.  At  the  present  time,  a  model  with  900,000  grid  points  can  be 
simulated  with  the  state-of-the-art  computer  system,  but  certain  tradeoffs 
have  to  be  made.  Simple  and  less-effident  numerical  integration  schemes 
have  to  be  employed  so  that  only  limited  neighboring  computational  points 
reside  within  the  computer’s  main  memory  at  a  given  time.  Extensive 
data-swapping  has  to  take  place  between  the  CPU  and  the  external  memory 
units.  As  a  tradeoff,  a  simulation  will  take  a  longer  time  depending  on  the 
numerical  scheme  and  the  computational  hardware  involved. 

In  a  CFD  model,  the  friction  coefficient  is  usually  expressed  as  a  function 
of  the  computed  velocity  profile  perpendicular  to  the  wall,  whereas  the  heat 
transfer  coefficient  is  usually  expressed  as  a  function  of  the  computed 
temperature  gradient  perpendicular  to  the  wall.  Consequently,  it  is  important 
to  provide  proper  boundary  conditions  at  the  wall  so  that  the  effect  of 
interdependency  is  minimized  during  the  computation  of  the  aerothermal  field 
around  the  NASP. 

For  the  past  35  years,  the  speed  and  memory  of  the  most  powerful  U.S. 
computers  have  increased  an  order  of  magnitude  every  seven  years.  The 
present  speed  limit  of  the  arithmetic  unit  is  between  one  and  two  billion 
floating  point  operations  per  second  (FLOPS)  depending  on  the  degree  of 
vectorization  of  the  code. 

At  present,  a  typical  computer  system  used  for  large-scale  CFD  simulation 
has  approximately  16  million  words  of  main  memory  with  possibly  one  or 
more  external  memory  units,  typically  8-mega  words  each.  Supercomputers 
are  nearly  all  vector  machines.  On  the  average,  the  relative  improvement  in 
speed  of  a  typical  CFD  code  is  approximately  six  times  after  vectorization. 

Beginning  in  the  spring  of  1989,  CFD  codes  are  expected  to  be  simulated 
on  CRAY’S  Y-MP/832,  which  offers  two  to  three  times  the  performance  of 
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the  X-MP.  The  system  has  8  CPUs,  operates  on  a  6-nanosecond  clock  cycle, 
and  has  a  central  memory  capacity  of  32  million  64-bit  words.  The  next  wave 
of  supercomputers  being  developed  by  companies  such  as  CRAY  and  NEC  will 
be  five  times  faster.  They  are  expected  to  be  introduced  within  five  years. 

Since  there  are  more  unknowns  than  just  the  number  of  (Navier  Stokes) 
equations  for  turbulent  flows,  it  is  necessary  for  CFD  models  to  "close"  the 
information  gap  using  turbulence  modeling.  Theoretically,  higher-order 
turbulence  models  are  more  universally  applicable  than  the  lower-order 
models  under  various  aerodynamic  conditions.  However,  this  has  not  always 
been  the  case  in  hypersonic  applications.  Most  pubbshed  models  for  NASP 
application  use  a  zeroth-order  closure  scheme  (see  Table  4,  Sec.  5,  for 
example).  This  particular  low-order  scheme  can  make  better  predictions  than 
many  higher-order  models  in  benchmark  tests.  Simpler  models  (usually  with 
fewer  constants),  if  designed  by  experienced  aerodynamicists,  can  sometime 
outperform  higher-order  models  with  more  "closure  constants."  Therefore, 
lab  and  design  experiences  also  play  an  important  part  in  CFD  modebng.  At 
present,  areas  with  low  predictive  reliability  in  turbulence  modeling  that 
should  be  emphasized  in  future  research  include 

1.  Improvement  in  prediction  of  the  location  and  length  of  the 
laminar-turbulent  boundary  transition  zone,  particularly  in  the  upper 
hypersonic  flight  range.  The  uncertainty  in  the  transition  process 
translates  directly  into  uncertainties  in  the  prediction  of  NASP 
performance,  weight,  and  other  control  parameters.  It  also  influences 
the  ability  to  predict  engine  inlet  characteristics. 

2.  Validation  of  higher-order  turbulent  closure  models  to  reduce  the  level 
of  uncertainties  in  the  area  of  low  predictive  reliabilities,  such  as: 

•  Strong  aerodynamic  curvature. 

•  Intermittency  and  large  structures  in  the  flow. 

•  Rapid  compression  or  expansion. 

•  Kinematically  influenced  chemical  reaction. 
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•  Low  Reynolds  number  effects. 

Strong  swirl. 

.  Turbulence,  which  is  strongly  influenced  by  body  force  acting 
in  a  preferred  direction. 

•  Uncertainties  in  boundary-condition  specification. 

•  Compressibility  effects  on  turbulence. 

■  Dynamic  stabibty  and  intermittency  in  the  high-wave¬ 
number  turbulent  eddies  in  the  SCRAMJET  combustion  process. 

•  Improvement  in  the  modebng  of  compressible  turbulence,  particularly 
near  the  engine  inlet  where  boundary  layer  and  shock  layer  may  induce 
substantial  density  variability  and  unstabibty. 

3.  Improvement  in  conserving  mass,  momentum,  and  energy  in  the 
numerical  schemes  for  solving  the  Navier  Stokes  equations  involving 
nonorthogonal  grid  transformation.  Need  effective  methods  for 
handling  sharp  gradient  and  discontinuity. 

4.  Improvement  in  treatment  of  boundary  conditions,  particularly  in  an 
interdependent,  nested  modeling  system  in  which  different  grids  are 
used.  Need  effective  numerical  scheme  for  handbng  nonreflective 
boundary  conditions  associated  with,  for  example,  a  three-dimensional 
full  Navier  Stokes  (FNS)  numerical  solver.  Reverse  flow  must  be 
addressed  effectively. 

5.  Development  of  quantitative  estimates  of  reductions  in  computational 
accuracy  near  the  nose  and  leading  edges  when  impUcit  schemes  (e.g., 
alternating  direction  implicit)  are  combined  with  the  coordinate 
tranformation. 

6.  In-depth  analysis  of  various  hypotheses  associated  with  turbulence 
modeling  and  the  universality  of  turbulence  closure  constants. 

7.  Development  of  methodologies  for  the  quantitative  assessment  of 
uncertainties  in  the  hypersonic  CFD  simulations.  It  is  important  that 


the  uncertainties  in  the  simulation  results  be  estimated  and  published  as 
a  function  of  the  vehicle  speed  and  location  along  the  vehicle  so  that 
designers  and  policymakers  can  make  reasonable  assessment  of  the 
vehicle’s  projected  performance. 

There  is  no  need  for  "turbulence  closure"  if  a  prototype  aircraft  is 
represented  by  a  model  containing  the  number  of  grid  points  approaching  the 
9/4^^  power  of  the  highest  flight  Reynolds  number.  For  the  anticipated  range 
of  flight  Reynolds  number  of  NASP,  say  10^,  it  would  need  a  computer  at 
least  10  times  more  powerful  than  the  pr^ent  system.  Consequently, 
turbulence  modeling  will  still  be  needed  for  the  foreseeable  future.  The 
present  trend  is  to  build  higher-order  "stre-'s-component"  type  models.  The 
simplest  stress-component  model  sol"  j^en  extra  partial  differential 
equations  in  addition  to  the  basic  bi^vier  Stokes  equations  of  motions  and 
continuity. 

Many  of  the  present  CFD  models  are  based  on  the  "parabobzed"  Navier 
Stokes  equations  (PNS),  which  makes  the  numerical  solution  process  much 
simpler  and  more  efficient.  The  parabolization  technique  neglects  the 
convective-diffusive  process  parallel  to  the  aircraft  surface.  For  NASP 
applications,  many  layers  are  generally  used  to  resolve  the  normal  gradient  by 
means  of  stretched  coordinates,  so  that  parabolization  is  often  justified. 
However,  comparison  of  PNS  and  FNS  solutions  should  be  made  to  evaluate 
the  importance  of  the  eliminated  convective  terms  in  the  governing  equations 
to  justify  the  use  of  a  PNS  code. 

The  majority  of  the  present  CFD  models  for  hypersonic  aircraft  use  the 
finite-difference  scheme  schematized  over  a  boundary-fitted  coordinate 
system.  Coordinate  transformation  becomes  one  of  the  major  efforts  in  CFD 
modeling.  Combination  of  the  (simpler)  finite-difference  method  (FDM)  and 
the  coordinate  transformation  replaced  the  finite-element  method  (FEM), 
which  enjoyed  some  popularity  during  the  1970s  because  of  its  geometric 
flexibility.  Since  the  finite-element  models  need  to  invert  extremely  large 
matrices  involving  the  variables  of  the  entire  aircraft,  they  are 
computationally  much  less  efficient  than  the  finite-difference  models.  The 
present  trend  is  to  use  the  most  straightforward  numerical  scheme,  so  that 


only  a  few  points  reside  within  the  CPU.  By  doing  so,  a  hypersonic  aircraft 
can  be  schematized  into  a  very  large  number  of  grid  points. 

In  terms  of  modeling  accuracy,  with  the  maximum  allowable  spatial 
resolution,  the  present  accuracy  of  CFD  models  at  the  lower  hypersonic  range 
(M  ^  6)  is  5  to  7  percent  in  terms  of  pressure.  When  the  predicted  values  are 
compared  with  the  experiments,  the  larger  deviations  are  located  near  the 
sharper  geometric  transitions.  This  indicates  that  higher  spatial  resolutions 
will  likely  improve  the  simulation  results.  A  higher-order  turbulence  model 
may  also  improve  the  results  around  the  curved  surfaces.  Both  imply  the  need 
for  bigger  and  faster  computers. 

As  far  as  hypersonic  laminar-turbulent  transition,  it  may  be  years  before 
we  have  either  the  data  or  the  model  to  go  beyond  crude  semi-empirical 
models.  Similarly,  not  enough  is  known  about  compressible  turbulence  to 
establish  suitable  closure  models  for  hypersonic  flows,  including  the  possibility 
of  large  turbulent  structures  or  perhaps  even  relaminization  from  a  turbulent 
state. 

Judging  by  past  trends,  by  the  year  2000  the  expected  speed  of  the  fastest 
computer  for  hypersonic  aircraft  applications  will  be  in  the  neighborhood  of 
10-50  billion  floating  point  operations  per  second.  However,  if  the  heat 
dissipation  problem  is  drastically  reduced  by  superconductive  material,  the 
potential  CPU  speed  may  substantially  improve  in  the  future. 

Because  of  the  lack  of  suitable  verification  data  at  high  speed,  CFD  codes’ 
ability  to  predict  depends  to  a  great  extent  on  the  universality  of  their 
turbulence  models.  In  other  words,  the  values  of  the  "closure  constants"  used 
in  these  models  would  have  to  stay  the  same  within  their  predictive  range. 
However,  several  well-known  universal  constants  in  fluid  mechnics  (e.g.,  the 
Kolmogorov’s  universal  constant)  were  found  to  be  variable.  In  Kolmogorov’s 
case,  when  a  3-D  turbulence  is  reduced  to  a  2-D  turbulence,  a  localization 
adjustment  factor  is  needed  when  computing  the  spectral  distribution  of 
turbulent  energy  near  a  wall.  Hypersonic  modeling  may  be  affected  if 
parabolization  is  involved. 

At  present,  major  hypersonic  CFD  research  is  conducted  in  government 
research  centers.  To  apply  CFD  methodologies  effectively  in  hypersonic 
aircraft  design,  these  modelers  not  only  have  to  work  with  experienced 
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designers  in  the  industry  but  also  have  to  work  with  material  engineers  to 
extend  fluid-dynamic  models  into  "aerothermoelastic"  design  codes. 
Furthermore,  it  is  urgent  that  an  effective  government /industry  relationship 
be  established  at  the  technical  level  so  that  research  codes  developed  in 
government  labs  can  become  production  code  that  designers  can  use.  It  is  also 
important  that  the  uncertainties  in  the  CFD  simulation  results  be  estimated 
and  published  as  a  function  of  the  vehicle  speed  and  location  along  the  vehicle 
(tip  to  tail),  so  that  the  NASP  design  teams  and  policymakers  can  estimate 
the  consequences  of  these  uncertainties  on  vehicle  performance. 
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Crj  c^  =  turbulence  modeling  constants 

CNS  =  compressible  Navier  Stokes 

Cj  =  tubulent  closure  constant  for  dissipation 

Cjj  =  drag  coefficent 

C^  =  Courant  number 

C|  =  skin  friction  coefficient 

CFD  =  computational  fluid  dynamics 

Cj^  =  lift  coefficient 

CPU  =  central  processing  unit 
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D 

DCa 

DEC 

del 


e 

ETAIO 


FDS 
PCX 
FFT 
F,  G, 

^Kleb 

FLOPS 

FNS 

FVS 

F 

wake 


^11, 

gy 

h 

«o 

BFF 

EGG 

HGV 


=  body  diameter 
=  binary  diffusion  coefficient 
=  dual  combustor  EAHJET 
=  Digital  Equipment  Corp.  (VAX  machines) 

=  deformation  operator 
=  diffusion  coefficients  in  a  coupled 

dynamic/chemistry  model 
=  total  SGS  energy  per  unit  mass 
=  the  latest  model  of  CYBER  series  supercomputer 
(ETA  Corp,  production  stopped  in  spring  1990) 

=  interfacial  friction  between  two  fluids 

=  flux  difference  splitting 
=  flux- corrected  transport 
=  fast  Fourier  transform 
H  =  vector  flux 

=  a  turbulence  modeling  parameter  (zero- equation 

model)  involving  Klebanoff  intermittency  factor 
=  floating-point  operations  per  second 
=  full  Navier  Stokes 
=  flux  vector  splitting 

=  a  parameter  in  turbulence  modeling  (zero- equation 

model)  involving  the  difference  between  the  max  and 
the  min  total  velocity 

g22=  components  of  the  contravariant  metric  tensor, 
or  gJ^ 

=  acceleration  due  to  gravity 

=  enthalpy 
=  total  enthalpy 

=  Ames  hypervelocity  free- flight  facility 
=  hyperbolic  grid  generator 
=  hypersonic  glide  vehicle 
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I 

'sp 

I/O 

J 

JFM 

k 


K 

K 

n 

L 

LACE 

LaRC 

LDV 

Le 


LHS 

M 

m,  H 
NAS 
NASP 
NDV 

n,  N 

ODE 

P 

PARIS 

PBR 

PDE 

Pe 


=  integer  index 

=  inter- fluid  source  term  in  a  two- fluid  model 
=  specific  impulse 

=  input /output 

=  Jacobian  of  a  transformation 
=  Journal  of  Fluid  Mechanics 

=  sub- grid- scale  turbulent  kinetic  erergy,  ergs/unit 
mass 

=  a  subscript 

r—  H 

=  Knudsen  number,  =  1.26  V  7 -jg- 

=  characteristic  length 
=  liquid  air  cycle  engine 
=  Langley  Research  Center 
=  laser  Doppler  velocimetry 

=  Lewis  number,  which  represents  the  relative  rate 
of  diffusion  of  mass  and  heat 
=  left-hand  side  (of  an  equation) 

=  liquid  oxygen 

=  Mach  number 
=  integer  index 

=  numerical  aerodynamic  simulator 
=  national  aerospace  plane 
=  NASP- derived  vehicle 
=  integer  index 
=  number  of  species 

=  ordinary  differential  equation 
=  static  pressure 

* 

=  commands  (used  together  with  C  language  for 
increasing  processing  speed  in  connection  machines 
=  Ames  pressurized  ballistic  range 
=  partial  differential  equation 
=  Peclet  number  =  Pr  x  Re  (Prandtl  No.  x 
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Reynolds  No.) 

=  control  functions  in  grid  generation 
PNS  =  parabolized  Navier  Stokes 

Pr  =  Prandtl  number,  which  is  the  ratio  of  rate  of 
diffusion  of  vorticity  and  of  heat 
Pr^  =  turbulent  Prandtl  number 

q  =  heat  flux 

R  =  volume  fraction  of  two- fluid  model 

Rg  =  free  stream  Reynolds  number 

RHS  =  right-hand  side  (of  an  equation) 

RHYFL  =  Rockdyne  Hypersonic  Facility  (a  shock  tube) 

RMS  =  root- mean- square 

RNS  =  reduced  Navier  Stokes 

R^  =  mass  diffusion  term  such  as  in  a  species  equation 

Rj^  =  Richardson  number,  a  measure  of  vertical  dynamic 
stability 

R^,  =  critical  Richardson  number 

cr 

S  =  source  term  within  fluid 

S^  =  Schmidt  number,  which  represents  the  relative 

rate  of  diffusion  of  vorticity  and  mass 
SGS  =  sub- grid- scale 

SSD  =  solid-state  storage  device 

t  =  time 

T  =  static  temperature 

te  =  trailing  edge 

[  ]  =  transposed  form  of  a  vector 

TIM  =  time  iterative  marching 

TKE  =  turbulent  kinetic  energy 

TINS  =  thin- layer  Navier  Stokes 

TVD  =  total  variation  diminishing 

0  =  dependent  variable 

UNS  =  unsteady  Navier  Stokes 
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upwind  PNS  solver 

vector  and  array  syntax  translator  (for  CFD  code 
vectorization) 

velocity  components  in  the  x,  y,  z,  direction 
velocity  fluctuation  components 

diffusion  velocity  components  for  air  species  in 
the  X,  y,  X  direction 

weighting  function  in  a  moving  adaptive  grid  system 

transposed  form  of  air  chemistry  source  term 

an  experimental  aerospace  plane 

the  y- distance  at  which  F  ^  occurs,  in  turbulence 
'  max 

modeling 

law- of- the- wall  distance  (i.e., 

Cartesian  gradient  vector  operator 
grid  spacing 

characteristic  length  of  viscous  layer;  boundary 
layer  thickness 
Kronecker  delta 

eddy  viscosity  or  isotropic  part  of  turbulent  energy 
dissipation  in  a  turbulence  model 
second  viscosity  coefficient,  or  conductivity 
diffusion  coefficient  within  one  fluid  of  a 
two- fluid  model 

first  viscosity  coefficient,  dynamic  viscosity 
turbulent  eddy  visicosity 

kinematic  viscosity 
vorticity 

air  chemistry  source  term 
density 

density  of  the  lighter  of  the  two  fluids  in  a 
two- fluid  model 
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average  density 

variation  of  density 
fluid  property 

gas  specific  heat  ratio  (1,4  for  perfect  gas) 
constants  in  the  transformed  space,  x  = 

y  =  yU,V,(),  z  =  z(^,>7,0 

turbulence  modeling  constants 

body  thickness 

shear  stress  components 
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1.  INTRODUCTION 

Advances  in  computational  fluid  dynamic  (CFD)  methods  and  the 
increasing  capabilities  of  computers  will  play  a  crucial  role  in  the  development 
of  the  National  Aerospace  Plane  (NASP).  Because  it  is  not  yet  practical  to 
conduct  continuous  flight  experiments  at  speeds  much  above  Mach  8, 
aerodynamic  behavior  and  the  potential  performance  of  the  NASP  beyond  this 
speed  have  to  be  predicted  using  the  CFD  methods  without  sufficient 
verification  data.  Therefore,  the  degree  of  certainty  or  the  confidence  limits 
associated  with  CFD  predictions  for  various  NASP  applications  are  at  this 
moment  difficult  to  quantify  (Mehta,  1990b).  This  is  particularly  true  for  the 
newly  developed  computational  codes  based  on  the  full  Navier  Stokes  (FNS) 
equations  with  finite-rate  chemistry.  It  is  a  challenging  task  to  prove  the 
validity  of  supersonic  combustion  processes  and  the  performance  parameters 
associated  with  the  SCRAMJET  engine  at  upper  hypersonic  speed  ranges 
using  a  computer  simulation. 

Since  CFD  models  are  based  on  numerical  solutions  of  a  set  of  governing 
equations  (under  certain  assumptions  over  a  finite-grid  network),  to  simulate 
the  aerospace  plane’s  integrated  external  flow  field  and  propulsion  system,  a 
numerical  model  must  cover  several  different  flight  regimes.  Although  much 
progress  has  been  made  in  solving  aerodynamic  design  problems,  many  new 
developments  are  still  needed  before  the  complete  three-dimensional  equations 
for  unsteady  compressible  viscous  flow  can  be  solved  routinely  with  high 
enough  resolution  and  certainty  in  the  entire  operating  range  of  the  NASP. 

OBJECTIVE 

The  research  reported  in  this  Note  sought  to  evaluate  independently  the 
degree  of  uncertainty  and  the  technical  development  risk  involved  in 
predicting  NASP  performance  using  numerical  simulation  of  the  highly 
integrated  air  frame/propulsion  system.  This  Note  presents  a  technical  review 
and  a  summary  of  four  interrelated  fields. 
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ORGANIZATION  OF  TfflS  REPORT 

Section  2  gives  the  hypersonic  environment  within  which  the  NASP  is 
expected  to  operate.  It  summarizes  governing  dynamic  equations  and  reviews 
many  forms  of  the  reduced  equation  sets.  Section  3  reviews  the  important 
numerical  schemes  upon  which  CFD  codes  are  based.  It  discusses  advantages 
and  disadvantages  as  weU  as  the  computational  requirement  associated  with 
each  numerical  scheme.  It  also  covers  one  of  the  important  characteristics  of 
the  hypersonic  aerodynamic  process  -the  existence  of  a  shock  wave  in  the  flow 
field-  as  well  as  the  use  of  a  numerical  technique  and  finite  computational  grid 
to  resolve  this  shock  wave. 

Section  4  addresses  a  new  technical  field  in  computational  fluid  dynamics: 
the  grid  generation  technique.  An  essential  as  well  as  a  time-consuming  step 
in  CFD  modeling  is  to  generate  a  proper  grid  around  an  air  frame  or  through  a 
duct  system.  Theoretically,  each  time  a  numerical  integration  step  is  carried 
out  over  a  transformed  grid  system,  a  certain  amount  of  mass  will  be  lost. 

The  aspects  of  conservation  of  mass,  energy,  and  momentum  associated  with  a 
transformed  computational  domain  are  also  discussed  in  this  section. 

Section  5  discusses  the  computational  needs  for  various  types  of  turbulence 
modeling.  Inasmuch  as  the  speed  and  the  size  of  present-day  computers 
cannot  yet  represent  the  flow  field  with  high  enough  resolution,  we  still  need 
to  model  sub-grid-scale  turbulence.  Turbulent  flow  will  occur  within  the 
flight  regime  of  the  aerospace  plane  and  the  supersonic  combustion  process  of 
the  SCRAMJET  engine.  The  modeling  of  turbulence  is  essential  in  predicting 
the  advective  process  of  momentum  and  heat  transfer  within  the  boundary 
layer  of  the  aerospace  plane.  This  in  turn  determines  the  potential 
performance  and  the  temperature  distribution  throughout  the  NASP  This 
section  also  describes  the  zeroth-,  first-,  and  second-order  as  well  as 
stress-component  turbulence  closure  schemes.  Computational  needs  for 
various  types  of  turbulence  modeling  are  discussed  in  this  section. 

Section  6  discusses  the  trends  in  development  of  computers  that  have  been 
used  to  carry  out  CFD  modeling.  To  use  a  CFD  model  effectively,  good 
hardware  and  software  support  are  essential.  Code  development  has  to  match 
the  type  of  machine  for  efficient  simulations.  Programming  and  CFD 
modeling  can  be  carried  out  more  effectively  by  using  proper  computers. 

Some  predictions  are  also  included. 
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Section  7  summarizes  CFD  modeling  needs  and  discusses  its  uncertainties. 
It  is  beyond  the  scope  of  this  Note  to  give  a  quantitative  analysis  of  technical 
risk  associated  with  predicting  the  NASP  performance  using  CFD  simulation. 
This  Note  does,  however,  point  out  several  aspects  pertaining  to  the  needs  in 
NASP  design  using  CFD. 

Section  8  presents  conclusions. 
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2.  HYPERSONIC  ENVIRONMENT  AND  THE  GOVERNING 
NAVIER  STOKES  EQUATIONS 


Numerical  models  designed  to  simulate  the  external  aerodynamics  or  the  duct 
flow  (of  SCRAMJET  for  examplej  of  an  aerospace  plane  should  rest  on  a  set  of 
appropriate  governing  equations.  These  equations  must: 

•  Include  essential  physical  parameters  in  the  basic  formulation, 

•  Operate  under  appropriate  assumptions  and  limits,  and 

•  Have  the  necessary  initial  and  boundary  conditions  for  the  computa¬ 
tion. 

To  make  a  reasonable  assessment  of  the  computational  needs,  the  range  of  the 
aerospace  plane’s  operational  parameters  is  discussed  first.  This  discussion 
includes  speed,  skin  temperature,  air-density  range,  and  other  physical  factors 
that  are  essential  to  the  plane’s  aerodynamics. 


THERMO-  AND  AERODYNAMIC  OPERATING  RANGE  OF  THE 
AEROSPACE  PLANE 

The  aerospace  plane  is  propelled  by  both  RAMJET  and  SCRAMJET 
air-breathing  engines  at  different  speed  ranges.  The  primary  difference  bevween 
the  two  engines  is  that  air  passes  through  the  RAMJET  at  subsonic  speeds  and 
through  the  SCRAMJET  at  supersonic  speeds.  In  a  RAMJET,  air  is  slowed 
down  to  about  Mach  0.2;  fuel  is  then  added.  At  hypersonic  speeds,  the 
SCRAMJET  is  designed  to  compress  and  decelerate  the  incoming  air  to  about 
Mach  1.  Adding  heat  to  the  duct  slows  air  further.  The  resulting  back  pressure 
creates  a  shock  wave  in  the  inner  diffusor  behind  the  inlet  through  but  ahead  of 
the  combustion  chamber.  At  around  Mach  6  (1800  m/s),  inlet  temperature  is 
approximately  2540  (Mackley,  1986).  At  that  speed,  the  inlet  air  speed  is 
slowed  to  about  1500  m/s.  Air  speed  at  the  hydrogen  injector  is  also  around  1500 
m/s.  Hydrogen  resides  in  the  combustor  for  approximately  one  milbsecond. 
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Unlike  conventional  aircraft,  the  fuselage  of  the  aerospace  plane  has  to  supply 
part  of  the  compressional  surface  for  the  engine.  The  same  is  true  for  the  aft 
portion  of  the  plane  and  exhaust  nozzle  for  propulsion.  Therefore,  the 
aerodynamic  computation  of  the  engine  and  fuselage  not  only  has  to  be 
dynamically  coupled  but  also  properly  posed,  otherwise  it  may  create  redundancy 
in  specifying  the  boundary  conditions. 

The  external  temperature  of  the  aerospace  plane  during  supersonic  cruise 
varies  at  different  flight  regimes  within  the  air-breathing  operational  range. 
Figure  1  superimposes  the  continuous  flight  corridor  (Masson  and  Gazley,  1956) 
over  the  aerospace  plane’s  flight  envelope.  The  maximum  ranges  of  external 
temperature  at  various  points  of  the  aerospace  plane  are  estimated 
(Korthals-Altes,  1987)  and  appear  in  Table  1. 

Table  1 

The  Maximum  Range  of  External  Temperature  at 
Various  Points  of  the  Aerospace  Plane 


_ Location _ 

Nose . 

Outside  of  the  crew  cabin . 

Fuselage  structure . 

Outside  passenger  and  cargo  area. 
Outside  liquid-hydrogen  tanks.... 

Wing  structure . 

Air  inlet  and  propulsion  unit . 

Control  surface . 

Leading  edges . 


Estimated  max,  temp. 

.  3260  °F 

....  2000  °F 

.  1800  °F 

....  1600  °F 

.  1400  °F 

....  2650  °F 

.  1600  °F 

.  3250  °F 

.  2650  °F 


The  equilibrium  skin  temperature  listed  here  is  the  temperature  at  which  the 
convective  heat  transfer  to  the  surface  equals  the  radiation  heat  transfer  from  the 
surface. 


0  5  10  15  20  25 

Mach  number 

SOURCE:  Adapted  from  Masson  and  Gazley  (1956). 

Fig.l — Assumed  Prototype  NASP/NDV  Characteristics  and  the  Flight 
Envelope  of  the  Aerospace  Vehicle 
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The  range  of  Reynolds  number  per  foot  of  characteristic  length  within  the 
flight  regime  is  from  10^  to  10®.  Boundary-layer  transition  from  laminar  flow  to 
turbulent  flow  will  probably  occur  along  the  hypersonic  body.  The  type  of  fluid 
flow  for  the  numerical  modeling  will  include  a  viscous  shock  layer  and  viscous 
boundary  layer.  In  the  air-breathing  flight  regime,  the  air-density  ratio  relative 
to  the  sea-level  value  and  the  mean  free  path  is  reflected  in  Table  2. 


Table  2 

The  Air-Density  Ratio  Relative  to  the  Sea-Level  Value 
in  the  Air-Breathing  Flight  Regime 


Altitude  (ft) 

Mean  Free  Path  (ft) 

Air  Density 

Relat  i  ve  1 0  Sea  Level 

0 

2.18  X  lO"”^ 

1 

100000 

1.61  X  10“^ 

1.35  X  10"^ 

200000 

8.45  X  10"^ 

2.57  X  10^ 

300000 

1.25  X  10"^ 

1.75  X  10"® 

From  Table  2,  the  Knudsen  number  within  the  aeroplane’s  air-breathing 

—7 

operating  range  is  between  2.18  X  10  to  around  0.1  per  foot  of  characteristic 
length  (body  dimension  or  boundary-layer  thickness).  In  a  strict  sense,  flow  with 
the  Knudsen  number  greater  than  unity  should  be  classified  as  rarefied  gas  flow 
which  is  not  in  a  continuum  regime.  Therefore,  it  might  not  be  properly 
described  by  the  compressible  Navier  Stokes  equations.  Certain  fully  merged 
layer  regimes  can  only  be  classified  as  near-continuum,  and  the  Navier  Stokes 
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equations  should  be  used  only  for  crude  estimates  of  the  local  skin  friction  and 
heat  transfer  over  that  area.  However,  experimental  evidence  indicates  that  the 
Navier  Stokes  equations  have  an  appreciably  wider  range  of  empirically  justified 
validity  than  their  theoretically  defensible  range.  Our  analysis  will  cover  the 
aerodynamics  of  the  air-breathing  flight  regime  of  the  aerospace  plane  that  can  be 
described  by  the  Navier  Stokes  equations. 


FLOW  COMPUTATION  IN  THE  SCRAMJET 

A  major  component  of  the  aerospace  plane  is  the  supersonic-combustion 
RAMJET  (SCRAMJET).  SCRAMJET  adds  heat  either  to  product  gases  or  by 
direct  injection  of  the  fuel  into  the  supersonic  air  stream  with  simultaneous 
mixing  and  burning  in  the  combustion  chamber  at  supersonic  flows.  The  ideal 
SCRAMJET  cycle  is  based  on  Rayleigh  flow  in  which  friction  losses  are  neglected. 
The  thermo-gas  dynamic  processes  involved  are  (a)  heat  addition  at  supersonic 
speeds  (Rayleigh  flow),  (b)  isentropic  nozzle  expansion  to  atmospheric  pressure, 
and  (c)  heat  rejection  at  constant  pressure  by  exhaust  to  the  atmosphere. 

A  unique  characteristic  of  the  SCRAMJET  engine  is  that  the  ratio  of  kinetic 
jet  energy  to  input  heat  ranges  from  approximately  1.5  to  3.0  (Ferri,  1964).  In 
fact,  a  SCRAMJET  is  the  only  propulsion  system  capable  of  generating  relative 
jet  kinetic  energy  to  input  energy  ratios  greater  than  unity  (turbojet’s  ratio  is  0.3 
and  rocket’s  0.4).  The  idealized  overall  efficiency  of  a  SCRAMJET  system  is 
generally  quite  high  and  increases  with  increasing  flight  Mach  numbers.  On  the 
other  hand,  phenomena  such  as  combustion-shock  wave  interaction  and  frictional 
effects  may  set  limitations  and  constraints  on  the  simplified  SCRAMJET  cycle. 
The  addition  of  heat  to  the  combustion  process  creates  combustion  waves 
(detonation  or  deflagnation).  Combustion  and  shock  wave  interaction  in  a 
supersonic  RAMJET  system  is  quite  complex,  particularly  if  frictional  effects 
and  chemistry  are  also  considered.  Most  important  is  the  question  of  whether 
adequate  mixing  can  occur  in  a  reasonable  length,  or  will  the  SCRAMJET  be 
inefficient  in  practice?  These  aspects  need  extensive  analysis  using  models 
containing  all  the  essential  variables. 

In  designing  a  SCRAMJET,  the  area  of  the  combustor  has  to  be  increased 
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slightly  to  prevent  thermal  choking.  But  the  amount  of  increase  must  be 
carefully  controlled  so  as  not  to  extinguish  the  flame.  The  overall  system 
efficiency  of  a  SCRAMJET  is  very  sensitive  to  internal  skin  friction  loss  (Waltrup 
et  al.,  1979;  Schetz  et  al.,  1987).  The  internal  flow  dynamics  of  a  SCRAMJET 
engine  can  be  characterized  by  a  strong  nonuniform  supersonic  inflow  followed  by 
a  nozzle  with  no  contraction  (throat)  as  in  rockets,  RAMJETs  or  jet  engines. 
Therefore,  to  simulate  the  flow  characteristics  within  the  SCRAMJET  accurately, 
a  numerical  model  must  be  able  to  resolve  the  strong  convective  components  with 
chemical  reaction  effects  present  in  significant  amounts  at  the  end  of  the 
combustor. 

Because  of  the  strong  convective  acceleration  and  possible  separation  within 
the  SCRAMJET  engine’s  combustion  system,  simplified  numerical  schemes  often 
used  for  NASP’s  external  flow  simulations  (e.g.,  PNS,  discussed  in  Sec.  3)  may 
not  be  suitable  for  SCRAMJET.  For  gas  chemistries,  a  numerical  model  should 
at  least  include  the  molecules  formed  with  H,  N,  C,  and  O.  The  potential 
performance  of  a  SCRAMJET  engine  depends  to  a  great  extent  on  the  calculation 
of  its  nonuniformity  of  flow,  and  the  turbulent  mixing/combustion  process.  These 
aspects  put  stringent  demands  on  the  turbulence  modeling  computational  resource 
and  require  sufficient  data  for  model  verification. 

Until  recently,  modeling  of  combustion  processes  in  a  SCRAMJET  engine  was 
often  assumed  to  be  mixing-controlled  combustion  (complete  reaction).  Based  on 
this  assumption  a  two-dimensional  unsteady  Navier  Stokes  code  (TWODLE)  has 
been  developed  at  NASA  Langley  for  simulating  SCRAMJET  combustion  flow  in 
the  NASA  modular  SCRAMJET  engine.  This  CFD  code  is  formulated  to 
compute  the  combustion  of  hydrogen/air  mixture  using  a  complete  reaction  model 
in  which  instantaneous  reaction  is  assumed  at  any  point  where  both  fuel  and  air 
are  present.  Boundary  conditions  at  the  wall  are  assumed  to  be  adiabatic  and 
noncatalytic  (White  et  al.,  1987).  Chemical  reactions  in  a  SCRAMJET 
combustion  process  are  generally  controlled  by  kinetics  and  not  mixing.  However, 
under  certain  conditions  it  may  be  that  both  reaction  time  and  mixing  time  are 
critical.  NASA  Langley,  together  with  Johns  Hopkins  University  Appbed  Physics 
Laboratory,  recently  developed  a  9-species,  18-equation  finite-rate  chemistry 
code  named  SPARK  to  simulate  the  hydrogen-air  combustion  to  describe  the 
important  reactions  of  various  SCRAMJET  designs  (White  et  al.,  1987). 
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Time-dependent  solutions  of  chemically  reacting  flow  in  SCRAMJET  are 
often  quite  "stiff  numerically.  In  other  words,  the  characteristic  time  scale 
associated  with  the  chemical  reactive  flow  varies  widely.  The  wide  range  in  time 
scale  comes  from  the  fact  that  to  resolve  the  fine  details  of  the  flow  field,  the 
computational  grid  needs  to  be  small  to  have  higher  spatial  resolution  thus 
resulting  in  a  very  small  time  step  (see  Sec.  3). 


GOVERNING  EQUATIONS  FOR  THE  THREBi-DIMENSIONAL  FLOWS 


The  three-dimensional  Navier  Stokes  equation  in  a  Cartesian  coordinate 
system  in  a  conservation-law  form  (Anderson  et  al.,  1984)  is  as  follows; 
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where  U  =  dependent  variable; 

F,  G,  H  =  flux  vector; 

E^  =  total  energy  per  unit  volume; 
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p  =  static  pressure; 
q  =  heat  flux; 

T  =  viscous  stress  components; 


The  complete  Navier  Stokes  equation  is  highly  nonlinear,  and  it  does  not 
possess  any  known  analytical  solution.  At  the  present  time  it  can  be  solved  only 
approximately  using  numerical  methods.  Solutions  to  the  Navier  Stokes  have 
many  applications  in  aerodynamic  design.  For  example,  the  skin  friction 
coefficients  for  a  given  design  can  be  computed  with  proper  initial  and  boundary 
conditions.  The  skin  friction  coefficient  Cj  can  be  defined  as: 


/  5u_s 
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Rg  denotes  the  free  stream  Reynolds  number.  The  derivative  represents  the 

velocity  profile  normal  to  the  surface.  These  velocities  are  the  computed 
velocities  resolvable  by  a  model’s  grid  network.  Similarly,  a  local  heat  transfer 
coefficient  can  also  be  estimated  as  a  function  of  the  computed  local  temperature 
gradient. 

Models  based  on  the  complete  set  of  Navier  Stokes  equations  require 
substantial  computational  resources.  To  improve  the  computational  efficiency,  the 
variables  on  the  continuous  physical  space  are  first  mapped  onto  a  discrete, 
transformed  network  for  modeling.  For  most  practical  applications  in  hypersonic 
aerodynamics,  a  thin-layer  approximation  further  improves  the  computational 
efficiency. 


PARABOLIZATION  AND  THE  REDUCED  FORMS  OF 
NAVIER  STOKES  EQUATIONS 

Even  though  the  governing  Navier  Stokes  equations  provide  an  excellent 
description  for  fluid  flow,  they  are  of  an  elliptic  type  for  which  the  numerical 
techniques  are  lengthy  and  cumbersome.  They  often  require  iterative  solutions. 
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The  presence  of  thin  shock  and  shear  layers  makes  it  difficult  and  costly  to  solve 
within  an  acceptable  level  of  accuracy.  But,  if  one  can  neglect  aU  the 
stream-wise  viscous  diffusion  terms  and  modify  the  stream-wise  convective  flux 
vector  to  permit  stable  time-like  marching  of  the  equations  downstream  from  the 
initial  data,  then  mathematically  the  nature  of  the  Navier  Stokes  equation 
changes  from  elliptic-type  to  parabolic-type.  The  equations  can  then  be  solved 
by  advancing  an  initial  "plane"  of  data  in  space  rather  than  by  advancing  an 
initial  "cube"  of  data  in  time.  The  modeling  efficiency  drastically  increases. 

From  the  physical  point  of  view,  parabolization  neglects  the  diffusion  process 
parallel  to  a  surface  but  retains  all  three  momentum  equations  and  makes  no 
assumption  about  the  pressure.  The  diffusion  term  involving  derivatives  parallel 
to  the  surface  for  the  high  Reynolds  number  turbulent  flows  near  the 
vorticity-generating  surface  is  usually  extremely  small.  Many  layers  are  needed 
to  resolve  the  normal  gradient  by  means  of  stretched  coordinates,  so  that  the 
thin-layer  approximation  is  often  justified.  In  fact,  many  of  the  hypersonic 
aerodynamic  computational  schemes  presently  available  are  of  the  parabolized 
Navier  Stokes  (PNS)  type.  Detailed  derivation  of  the  PNS  equations  from  the 
complete  Navier  Stokes  equations  can  be  found  in  Rudman  and  Rubin  (1968)  or 
in  text  books  such  as  Anderson  et  al.  (1984). 

There  are  several  other  simpUfied  forms  which  deviate  slightly  from  PNS. 
Barnett  and  Davis  (1985)  added  the  stream-wise  dependence  to  PNS  by  means  of 
a  pressure  relaxation  scheme.  This  reduced  Navier  Stokes  (RNS)  method  is  valid 
for  air  flows  in  which  the  stream-wise  interaction  is  weak.  The  same  is  true  for 
the  thin  layer  Navier  Stokes  equations  (TLNS)  introduced  by  Pulliam  (1984).  In 
TLNS,  the  viscous  terms  containing  derivatives  in  the  directions  parallel  to  the 
body  surface  are  neglected  in  the  unsteady  NS  equations,  but  all  other  terms  in 
the  momentum  equations  are  retained.  By  retaining  the  terms  that  are  normally 
neglected  in  boundary^ayer  theory,  separated  and  reverse  flow  regions  can  still 
be  calculated  (Anderson  et  al.,  1984). 

Numerical  models  based  on  the  simplified  Navier  Stokes  equation  have  played 
a  key  role  in  the  design  of  future  aerospace  planes  because  they  reduce  the 
computational  effort  by  at  least  an  order  of  magnitude.  On  the  other  hand,  a 
simplification  scheme  such  as  the  "zonal  modeling"  technique  may  also  create 
potential  computational  problems.  For  example,  there  are  potential 
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incompatibilities  in  establishing  the  mutual  boundary  conditions  for  the  patched 
models,  particularly  when  the  flow  fields  of  these  models  are  dynamically  coupled, 
as  in  the  case  of  the  aerospace  plane. 

L'ke  many  other  types  of  the  Navier  Stokes-based  numerical  models, 
sometimes  substantial  efforts  are  required  to  establish  a  proper  boundary 
condition  so  the  interior  solution  is  not  influenced  by  the  reflective  boundaries. 
Theoretically,  computers  will  never  be  large  enough  that  the  model  boundaries 
can  be  placed  to  represent  the  truly  far  field  unless  we  know  the  solution 
beforehand.  To  establish  a  dynamically  suitable  non-reflective  boundary  in  a 
moving  coordinate  system  with  shock  is  understandably  quite  complex. 

Past  experiences  m  Navier  Stokes  modeling  have  indicated  that  more  accurate 
numerical  sc.  .es  generally  need  better  treatment  at  the  boundaries.  This  is 
because  ’jf  ur-order  (more  accurate)  numerical  schemes  involve  higher  spatial 
derivatives  that  need  to  be  specified  at  the  boundaries  during  the  integration 
process.  As  a  consequence,  models  based  on  different  solution  schemes  require 
different  boundary  treatments.  They  will  be  discussed  together  with  the  schemes 
themselves  in  the  Sec.  3. 

THREE-DIMENSIONAL  SYSTEM  IN  A  TRANSFORMED  SPACE 

The  shape  of  most  aerodynamic  bodies  consists  of  curved  surfaces. 
Consequently,  the  Navier  Stokes  equation  formulated  on  a  Cartesian  coordinate 
system  loses  its  resolving  power  in  aerodynamic  applications.  To  improve 
computational  efficiency,  it  is  more  desirable  to  formulate  the  equation  in 
curvilinear  coordinates.  An  arbitrary  shape  in  the  physical  domain  can  be 
mapped  onto  a  rectangular,  discrete  computational  domain  through  coordinate 
transformation.  The  transformed  space  is  usually,  but  not  necessarily,  a 
rectangle.  It  could  be  a  circular  disk  or  a  sphere  such  as  the  earth  ellipsoidal 
coordinate  (Liu  and  Leendertse,  1987).  The  great  majority  of  the  coordinate 
system  generation  methods  are  based  on  solving  partial  differential  equations 
using  a  coordinate  transformation  from  the  physical  domain  to  an  orthogonal 
computational  domain  (see  Fig.  2).  Using  this  approach,  local  refinements  of  the 
mesh  can  be  achieved  by  introducing  an  appropriate  mesh  control  function  into 
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the  governing  equations.  Boundary  fitted  coordinate  generation  methods  are 
discussed  in  Sec.  4. 

Three-dimensional,  mass-averaged,  Navier  Stokes  equations  in  the 
transformed  space  (^,ry,0  for  the  dependent  variable  U,  are: 

rar.  rdF.  rdF. 
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where  an  eddy  viscosity  model  is  assumed  and  where: 

=  Kronecker  delta; 

del  =  deformation  operator; 
f  =  eddy  viscosity; 

Pr  =  Prandlt  number; 

Pr^  =  turbulent  Prandtl  number; 

F,G,H  =  vector  flux  ; 

U  =  dependent  variables; 
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^,T],C  =  constants  in  the  ((,tiX)  space,  x=x(^,7/,0,  y=y{(,‘n,0>  etc. 

For  the  orthogonal  three-dimensional  transformed  system,  g^2  =  613  =  §23  ’ 

gjj  =  0  (i^j). 

The  quantity  g^j  in  curvilinear  coordinates  specifies  the  geometry  considered  and 

constitutes  the  matrix  tensor  of  that  space,  g-j  is  proportional  to  the  cosine  of  the 

angle  between  a  coordinate  line  along  which  the  curvilinear  coordinate  ^  varies 
(the  other  two  coodinates  being  constant  along  such  a  line). 

The  Laplacians 

=  (Siig22633)  ('^822^33^^11^ 

=  (811822833^  ^  ('^811633/622^ 
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It  must  satisfy  the  condition  within  the  flow  field: 
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or 
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Two-dimer.sional  orthogonal  systems  must  satisfy  these  conditions  on  the 
boundaries: 


-  ’!yVg22/6ll  iy  -  ’!x'^S22/®U 

=  -y^Vg22/gii  y,  =  *f''g22/gn 


If  the  system  is  orthogonal,  it  has  the  least  extra  terms  in  the  transformed 
PDE  to  account  for  the  curvature  of  transformation  and  the  centrifugal  force  of 
the  flow.  One  characteristic  of  the  orthogonal  system  is  the  vanishing  of  the 
off-diagonal  elements  of  the  matrix  tensor. 

Therefore,  the  basic  orthogonally  transformed  equations  are,  if  we  define 
Jacobian  of  the  coordinate  transformation, 
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which  measures  the  area  of  a  cell  (volume  in  a  three-dimensional  system), 
v/gji/gjj  =  cell  aspect  ratio  which  measures  the  ratio  of  the  length  of  the  sides . 


THREE-DIMENSIONAL  DYNAMIC  SYSTEM  COUPLED  WITH 
NONEQUILIBRIUM  CHEMISTRY 


When  the  aerospace  plane  reaches  high  altitudes  where  the  atmosphere  is 
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characterized  by  low-density  effects  and  reacting  gas  chemistry,  to  simulate  the 
flow  field  properly  the  governing  reacting  gas  chemistry  equations  have  to  be 
dynamically  coupled  to  the  Navier  Stokes  equations.  Recently  Hoffman  et  al. 
(1988)  made  some  tests  using  the  simple  Baldwin/Lomax  zero  equation 
turbulence  model  (Sec.  5)  in  the  full  Navier  Stokes  equation  formulation 
dynamically  coupled  with  nine  chemical  equations  for  air.  According  to  Hoffman 
et  al.,  nine  species  (include  electron  density)  have  been  found  to  be  important  in 
the  dynamic  coupling.  However,  when  the  equation  set  is  solved  using  an  implicit 
numerical  scheme  to  remove  the  chemistry  stiffness  problem,  the  size  of  the  block 
matrix  becomes  very  large.  For  example,  even  with  the  simplest  turbulent 
closure  scheme  (i.e.,  zeroth-order  scheme,  see  Sec.  5),  the  size  of  the  block  matrix 
becomes  14  x  14  in  which  five  are  Navier  Stokes  components.  Using  the  strong 
conservation  form,  in  Cartesian  space  the  equation  set  is  as  follows: 

d  V’  d  F’  d  G'  d  El 

- - C_^  yj 

d  t  d  X  d  y  d  z 

where 

T 

Up  =  . ' 

Pp  /?2> . ’  ~  density  of  the  species 

Ng  =  number  of  species 

Ns 

p  =  total  density  =  E  n 

s=l  ^ 

u,  V,  w  =  Cartesian  velocity  components 

rp 

W  =  chemistry  source  term,  =  [wp.,.  ,  ,  0,0, 0,0] 

s 

£4;p..Wj,j  =  source  term  for  each  species 

F^,  G^,  H^  =  coupled  dynamics/chemistry  flux  vectors 
each  with  an  invicid  and  a  viscous  part. 
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In  each  of  the  above  vectors  F’^,  G’^,  ,  the  left-hand-side  member 

represents  the  invicid  part  and  the  right-hand-side  member  represents  the 
viscous  terms. 

Furthermore,  Uj,...,  Uj^  are  diffusion  velocities  for  species  1  through  in 
^  s 

the  X  direction.  Similarly,  v  and  w  are  diffusion  velocity  components  in  the  y  and 

z  direction.  In  numerical  simiilations,  these  diffusion  velocity  components  are 

often  assumed  to  be  proportional  to  the  local  concentration  gradient  of  the 

species.  To  estimate  the  rate  of  diffusion,  diffusion  coefficients  such  as  D  are 

0 

used. 

In  a  strict  sense,  to  model  chemical  reactions,  one  has  to  consider  the 
nonequilibrium  internal  energy  modes.  The  total  internal  energy  is  the  sum  of 
the  energy  due  to  translation,  rotation,  vibration,  chemical 
disassociation/association,  electronic  state,  and  the  kinetic  energy  in  the  air  flow. 
However,  for  the  Navier  Stokes  equation  to  be  valid,  it  is  often  assumed  that  the 
translation  and  rotational  modes  are  in  equilibrium. 

A  significant  portion  of  the  energy  of  air  behind  the  bow  shock  may  be 
associated  with  the  dissociation  of  molecules  into  atoms.  Dissociation  effects 
begin  at  about  Mach  10.  When  the  chemical  reaction  rates  are  fast  (i.e., 
equilibrium),  the  heat  transfer  process  is  independent  of  the  mode  of  transfer  (by 
conduction  or  diffusion).  However,  in  modeling  high-speed,  low-density  flows, 
the  boundary  conditions  at  the  wall  often  have  to  consider  the  efficiency  of  the 
wall  in  catalyzing  the  recombination  of  atoms  into  molecules.  The  efficiency  has 
a  significant  effect  on  heat  transfer.  Chemical  kinetics  of  catalytic  surface 
reactions  in  a  hypersonic  viscous  shock  layer  of  a  nonequibbrium  gas  has  been 
studied  by  Chung  (1961).  Wray  (1962)  considered  seven-species  ionization  and 
Bittker  and  ScuUin  (1984)  and  Hoffman  et  al.  (1988)  formulated  a  nine-species 
air  chemistry  reaction  model  with  wall  catalysis. 

One  problem  of  modeling  low-density  real  gas  flow  is  the  "curse  of 
dimensionality,"  i.e.,  the  size  of  the  matrix  grows  rapidly  with  the  number  of 
species  being  considered.  This  is  particularly  true  when  higher-order  turbulence 
closure  schemes  are  used  in  the  model.  These  aspects  are  discussed  in  Sec.  5. 
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3.  NUMERICAL  SCHEMES 


Perhaps  the  most  important  factor  that  determines  the  computational 
requirement  in  modeling  is  the  principal  scheme  used  for  the  integration  of  the 
governing  Navier  Stokes  equation.  For  a  given  physical  problem  the 
numerical  scheme  determines  not  only  the  solution  accuracy  but  also  the 
requirements  in  computer  memory,  the  amount  of  arithmetic  operations,  and 
the  treatment  at  the  model’s  boundaries.  In  this  section  we  will  cover  the 
state-of-the-art  scheme,  the  most  popular  scheme,  and  those  with  the  highest 
potential  for  future  aerospace  plane  appUcations  considering  the  trend  of 
computer  developments.  Only  the  basic  numerical  scheme  will  be  treated 
here.  The  method  of  coordinate  transformation  is  discussed  in  Sec.4 


FINITE  DIFFERENCE  METHODS 
Explicit  Scheme 

The  original  explicit  finite  difference  scheme  proposed  by  MacCormack 
(1969)  was  the  primary  modeling  algorithm  for  nearly  a  decade  until  other 
relatively  more  efficient  schemes  became  more  popular.  However,  because  of 
its  simplicity,  the  basic  scheme  enjoys  extensive  use.  For  example,  the 
majority  of  recent  numerical  studies  involving  SCRAMJET  simulations  still 
use  this  simple  and  straight  forward  algorithm.  The  major  advantage  of  this 
method  is  its  ease  to  program  and  to  debug.  Its  major  disadvantage  is  the 
stability  conditions  that  limit  the  permissible  size  of  the  forward  integration 
time  step. 

The  MacCormack  explicit  scheme  is  a  predictor-corrector  type  of 
algorithm.  It  uses  a  split-operator  technique,  which  is  second-order  accurate 
in  time  and  space.  When  it  was  first  used  for  the  Euler  equation,  artificial 
diffusion  was  employed  thus  permitting  the  prediction  of  shocked  flows. 
Recently,  the  scheme  has  been  used  for  the  following  NASP  applications: 


Hypersonic  turbulent  mixing  and  reaction  of  hydrogen  fuel  and  air  near 

the  injectors  of  a  SCRAMJET  (Weidner  and  Drummond,  1982); 

The  SCRAMJET  flow  field  over  a  rearward  facing  step  with  a  transverse 

H  injector  that  includes  the  detail  binary  diffusion  of  hydrogen  and  air 
2 

along  with  variable  Lewis  number  (Berman  and  Anderson,  1983); 

The  SCRAMJET  flow  field  with  no  strut,  one  strut,  and  multiple  struts 
(Kumar,  1982); 

Wing  fuselage  aerodynamic  interaction  (Shang,  1984). 


The  explicit  numerical  algorithm  contains  essentially  a  predictor  and  a 
corrector  part.  For  illustration  purposes,  we  only  use  the  simplified  governing 
equation  (see  Eq.  (2.10))  in  the  conservation-law  form: 


mj  ,  dG 


=  0 


(3.1) 


The  u^^^  becomes  the  "current  value"  of  u  for  the  net  prediction  step,  and 
*  1 J 

so  on.  The  stability  conditions  associated  with  this  explicit  scheme  in  the  x 
and  y  direction  are,  respectively; 
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where 

c  =  speed  of  sound; 

7  =  ratio  of  specific  heat  of  gas; 

A,/i  =  viscosity  coefficients; 

Pr  =  Prandtl  number. 

For  NASP  applications,  a  typical  SCRAMJET  model  (e.g.,  30  x  109 
nodes)  using  an  explicit  scheme  takes  30,000  time  steps  to  reach  a 
steady-state  internal  flow  condition.  Using  a  first-generation  vector  machine 
such  as  the  Control  Data  Cyber  203  (equivalent  to  a  CDC  7600  non-vector 
machine),  it  takes  approximately  one  to  two  hours  for  the  numerical 
integration  process  (Weidner  and  Drummond,  1982). 

Vectorization  of  the  computer  code  can  also  improve  efficiency.  When  the 
explicit  scheme  is  vectorized  and  run  on  a  CRAY-1  vector  mechine  (rated  at 
160  million  floating  point  operations  per  second,  FLOPS),  the  vectorized 
program  outperforms  the  original  scalar  code  by  a  factor  of  8.31  (Shang, 

1984).  Using  the  original  explicit  scheme,  a  typical  wing-fuselage 
aerodynamic  interaction  simulation  (56,730  grid  points,  M=6)  requires  about 
1.4  hour  CPU  time  on  a  CRAY-1  (Shang,  1984). 

The  basic  explicit  scheme  has  also  been  modified  and  coupled  to  a  general 
interpolants  algorithm  (GIM,  Prozen  et  al.,  1977)  designed  to  take  advantage 
of  a  "pipeline"  feature  on  certain  computers  such  as  the  CDC  STAR  100 
machine.  To  compute  hypersonic  flows  (M  >  5),  a  factor  of  6  (in  speed)  was 
obtained  as  compared  to  the  equivalent  scalar  machine  (CDC  7600). 

In  all  of  the  above  computations;  the  explicit  scheme  is  coupled  to  the 
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algebraic  turbulence  models  proposed  by  Baldwin  and  Lomax  (1978)  and 
discussed  in  Sec.5 

Even  though  the  scheme  is  unconditionally  "stable"  according  to  linear 
stability  analyses  (Beam  and  Warming,  1976;  Pulliam  and  Steger,  1980), 
experience  indicates  that  the  scheme  with  centered  spatial  differences  has  only 
limited  range  of  Courant  numbers  (Thomas,  1988).  In  a  Cartesian 

context,  the  limiting  Courant  number  is 


At 


where 


for  perfect  gas 
where 

c  =  the  local  speed  of  sound 
Cq  =  the  Courant  number 

At  =  the  maximum  allowable  time  step 
u  =  the  vehicular  velocity 

I A I  =  the  largest  eigenvalue  of  the  invicid  flux  vector  Jacobian 
matrix. 

In  the  stream-wise  direction,  the  alternating  direction  implicit  (ADI) 
scheme  is  limited  to  a  Courant  nxunber  of  unity  (i.e.,  one).  In  the  direction 
normal  to  the  body,  a  much  higher  Courant  number  can  be  used.  In  addition, 
the  implicit  ADI  scheme’s  stability  is  sensitive  both  to  spatial  resolution  of 
the  grid  and  to  the  degree  to  which  the  instantaneous  solution  departs  from 
the  real  one  (Thomas,  1988).  To  satisfy  this,  the  initial  time  steps  for  the 
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integration  must  be  small.  Tassa  and  Conti  (1987)  suggested  using  different 
time  steps  for  each  grid  point  as  time  proceeds.  For  example,  for  a 
two-dimensional  flow,  to  determine  the  variable  time  step  size,  the  averaged 
eigenvalues  in  two  directions  have  to  be  evaluated. 


Implicdt-Elxplicit  and  Hybrid  Schemes 

When  used  for  viscous  or  turbulent  flows,  the  expbcit  method  described 
above  is  penalized  by  the  stiffness  of  the  discrete  Navier  Stokes 
approximation.  Sometimes  it  is  very  inefficient  in  simulating  flow  in  the 
flight  Reynolds  number  range.  In  the  explicit  method,  unknown  variables  are 
solved  using  known  variables  that  are  explicitly  defined  from  the  earlier  steps 
in  a  marching  manner.  In  the  implicit  methods,  however,  unknowns  are  first 
grouped  together  before  they  are  solved  simultaneously,  usually  by  means  of  a 
matrix  inversion  technique. 

Implicit  methods  can  often  be  designed  so  they  are  not  subject  to  certain 
linear  stability  conditions.  Therefore,  they  are  often  more  efficient  than  their 
explicit  counterpart  by  being  able  to  use  larger  permissible  time  steps. 

In  fact,  the  efficiency  of  aerodynamic  modeling  is  also  governed  by  other 
factors.  For  example,  at  high  Reynolds  numbers  the  magnitude  of  the  inertial 
force  described  by  the  hyperbolic  terms  of  the  Navier  Stokes  equations  is 
much  larger  than  the  viscous  force  described  by  parabobc  terms.  A  numerical 
scheme  can  be  made  more  efficient  by  splitting  the  equations  into  a  hyperbobc 
and  parabolic  part.  MacCormack  (1976)  proposed  a  method  that  solves  the 
hyperbolic  part  explicitly  by  using  characteristic  theory  and  solves  the 
parabolic  part  by  using  implicit  parabolic  method.  Both  methods  are  fully 
conservative  and  stable.  According  to  MacCormack,  under  the  thin  layer 
assumption,  the  velocity  component  normal  to  the  solid  surface  is  very  small 
compared  to  the  streamline  direction.  As  a  consequence,  the  stability 
condition  involving  the  normal  velocity  is  much  less  restrictive  than  the 
parabolic  part.  The  stability  condition  for  the  explicit  part  is: 
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<  -f4f-  (3-6) 

The  parabolic  part  solves  the  simple  tridiagonal  matrix  and  is 
unconditionally  stable.  The  semi-implicit  scheme  requires  roughly  1/6  of  the 
time  of  the  explicit  scheme.  Knight  (1984)  also  developed  a  hybrid  method 
that  combines  the  explicit  scheme,  described  earlier,  and  an  implicit  scheme 
for  the  viscous  sublayer  and  transition  wall  regions  of  the  turbulent  boundary 
layers.  When  the  hybrid  algorithm  is  coded  in  the  SL/1  vector  programming 
language  (developed  at  NASA  Langley),  the  computational  efficiency 
improved  by  a  factor  of  16  to  21  as  compared  to  a  vectorized  version  of  the 
MacCormack  time-split  alogorithm. 


Alternating-Direction  Implicit  Scheme 

The  most  adapted  multi-dimensional  implicit  method  is  the  ADI 
algorithm  introduced  by  Douglas  (1955),  Peaceman  and  Rachford  (1955),  and 
Douglas  and  Gunn  (1964).  The  ADI  scheme  was  modified  and  applied  for  the 
compressible  Navier  Stokes  equation  (Beam  and  Warming,  1976).  The 
scheme  is  second-order  accurate,  non-iterative,  and  linearly  unconditionally 
stable.  Even  though  the  scheme  is  a  three-time-level  scheme,  it  requires  only 
two-time-level  of  data  storage.  The  algorithm  allows  the  spatial 
cross-derivative  terms  to  be  included  efficiently  in  a  spatially  factored 
manner. 

Being  perhaps  the  most  popular  implicit  scheme  for  any  application, 
particularly  in  modeling  multi-dimensional  fluids,  it  consists  of  essentially  an 
explicit  and  an  implicit  portion.  Starting  from  the  specified  initial  field,  the 
two  portions  are  applied  in  an  alternating  manner  in  the  process  of 
integration. 

Several  codes  other  than  the  one  posed  by  Beam  and  Warming  (1976)  are 
based  on  the  general  ADI  scheme  (Steger,  1977;  Pulliam  and  Steger,  1980; 
Schiff  and  Steger,  1980).  The  generalized  ADI  scheme,  when  used  for  solving 
a  hyperbolic  nonlinear  two-dimensional  system  in  conservation-law  form,  is 
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where 

u  =  unknown  p-component  vector; 

F,  G  =  given  vector  functions  of  the  components  u; 
u(t)  =  u(nAt)  =  u°. 

i9F'  fiO 

Let  Jucobiun  matrices  A= — —  ;  F=  Au  und  G^Bu.  (3.8) 

The  second-order  recursion  algorithm  can  be  grouped  in  the  following 
sequences: 


-®j,k-l''j,k-lj 


(3.9) 


(3.12) 
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Predictor-Corrector  Type  Implicit  Scheme 

In  the  solution  process,  an  ADI  implicit  scheme,  as  discussed  in  Sec.  3, 
usually  involves  the  inversion  of  tri-diagonal  matrices.  Because  of  the 
recursion  and  inversion  procedures  associated  with  an  implicit  scheme,  it  may 
require  nearly  twice  as  many  operations  per  integration  step  as  an  explicit 
scheme.  Since  the  permissible  time  step  size  for  the  implicit  scheme  is  usually 
much  larger  than  the  explicit  scheme,  it  thus  justifies  the  increase  in  the 
number  of  operations.  It  would  be,  on  the  other  hand,  more  desirable  to 
design  an  implicit  scheme  that  reduces  the  number  of  operations  per 
integration  step.  MacCormack  (1982)  proposed  a  scheme  that  is  basically 
implicit,  but  it  allows  explicit  operation  in  the  region  away  from  the  aircraft 
where  mesh  spacing  is  larger.  The  method  is  based  on  a  predictor-corrector 
principle.  In  the  process,  it  solves  block  bi-diagonal  matrix  systems  rather 
than  the  block  tri-diagonal  one  associated  with  the  standard  ADI  scheme. 
Therefore,  for  the  overall  efficiency  for  hypersonic  applications,  this  scheme 
outperformed  the  pure  ADI  code  discussed  previously.  The  method  is  stable 
for  unbounded  At  and  second-order  accurate  under  the  condition  that  the 
following  term  remains  bounded: 


At  1 
Min(Ax^,At^) 


(3.13) 


5  7 

In  practical  computations  in  the  Re  range  between  3x10  to  3x10  , 
simulation  codes  based  on  this  method  can  make  stable  simulations  with 
Courant  number  as  high  as  1000  (MacCormack,  1982).  Considering  the 
governing  Navier  Stokes  equation  in  the  conservation  form  as: 


aj 

'Ux 


+ 


(3.14) 


The  equation  may  be  numerically  integrated  in  time  by  a 
predictor-corrector  sequence  outlined  as  follows: 
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1.  Prediction  step; 


AU“  .  =  -At 
1  .J 
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I-At— 
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2.  Correction  step: 
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(3.16) 


where  A  and  B  are  the  Jacobians  of  vectors  F  and  G  defined  by  A=5F /dU, 
B=dG/dU.  I  is  identity  matrix.  The  upstream  and  downstream  differencing 
operators  are  defined  by: 


A^Z.  .  Z. ,  ,  .  -  Z.  . 
+  _i,J_  1+1,  i  i,j 

Ax  Ax 


A  Z.  . 
-^= 


Z.  .  - 

__y_ 


7^ 


(3.17) 


for  variables  in  the  x  direction;  the  same  is  true  for  the  y  direction. 

Here,  only  the  essence  of  the  numerical  scheme  is  outlined.  Since  this 
scheme  may  be  the  highest  potential  for  the  aerospace  plane  application  in  the 
near  future,  several  of  its  operational  aspects  will  be  discussed  herein.  These 
include  stability,  memory  requirement,  and  the  relative  speed  as  compared 
with  the  other  schemes,  particularly  with  respect  to  the  traditional  ADI 
implicit  scheme  discussed  in  Sec.  3 
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Stability  limit:  Even  though  the  scheme  is  linearly  uncoTi<^itionaliy 
stable,  the  Navier  Stokes  equations  are  highly  nonlinear,  so  the 
stability  limit  for  various  applications  varies.  Numerical  experiments 
using  parabolized  codes  indicate  that  the  practical  limits  are  between 
Courant  number  of  100  to  1000.  The  maximum  limit  for  the 
predictor-corrector  explicit  scheme  is  around  1.5. 

When  the  Courant  number  exceeds  5000,  the  skin-friction  coefficients 
computed  using  this  method  become  inaccurate.  However,  on  fluid 
systems  with  free  surfaces,  the  order  of  accuracy  of  the  computed 
velocity  distribution  begins  to  drop  when  the  Courant  number  exceeds  3 
using  the  ADI  scheme. 

Memory  requirement:  The  predictor-corrector  scheme  requires  30 
percent  less  computer  memory  than  the  ADI  scheme  because  it  does 
not  need  to  store  the  entire  block  tri-diagonal  matrix  system.  Instead, 
it  forms  block  bi-diagonal  systems  and  inverts  with  one  sweep  while 
storing  only  two  4x4  matrices  at  a  time.  For  the  same  reason,  ADI 
schemes  involve  longer  computer  codes  than  this  scheme. 

Integration  speed:  Theoretically,  since  this  scheme  involves  only 
bi-diagonal  matrices  as  opposed  to  tri-diagonal  matrices,  it  should  be 
much  faster  than  the  ADI  scheme.  But  in  practice,  the  speed  gain  is 
somewhat  less  dramatic,  particularly  when  the  computational 
grid-network  is  evenly  spaced  so  that  only  the  implicit  mode  is  involved 
during  the  integration  process.  Another  major  factor  that  reduces  the 
speed  is  that  for  e/ery  integration  time  step,  the  |B|  matrix  on  the 
right-hand-side  (RHS)  has  to  be  inverted  twice.  The  same  is  true  for 
all  of  the  RHS  terms.  Strict  theoretical  analysis  indicates  that  two 
block  bi-diagonal  systems  can  be  inverted  about  10  percent  faster  than 
one  block  tri-diagonal  systems.  It  may  be  concluded  that  for  the 
predictor-corrector  type  of  implicit  scheme  to  be  more  efficient  than  the 
ADI  scheme,  the  former  must  contain  more  points  that  can  be 
computed  explicitly.  In  other  words,  the  grid  spacing  in  the  far  field 
should  be  much  larger  than  the  spacing  near  the  aircraft.  This  is 
usually  the  case,  however. 
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FINITB-ELEMENT  METHODS 

Before  the  finite-element  method  was  used  for  aerodynamic  computation, 
it  was  used  mainly  to  analyze  static  structure  systems.  As  early  as  1956  the 
method  was  employed  in  aircraft  frame  analysis.  Geometric  flexibility  has 
been  the  primary  feature  of  the  finite-element  scheme.  The  finite-element 
method  is  often  used  because  it  yields  good  accuracy  on  a  coarse 
computational  grid.  The  method  has  received  somewhat  less  attention  lately 
because  good  spatial  resolution  can  be  obtained  by  using  an  efficient  finite 
difference  scheme  coupled  with  a  coordinate  transformation  technique.  For 
NASP  applications,  the  finite-element  method  competes  directly  with  the 
finite  difference/coordinate  transformation  technique  as  the  primary 
aerodynamic  design  tool  in  the  higher  speed  range. 

The  finite-element  method  was  originally  formulated  for  linear  problems 
until  it  combined  the  method  of  weighted  residuals  (Finlayson,  1972)  and  the 
Galerkin  criterion  (Galerkin,  1915),  and  it  then  became  capable  of  handling 
problems  described  by  the  complete  nonlinear  partial  differential  equations 
such  as  the  Navier  Stokes  equations.  For  time  integration,  explicit  schemes 
always  seem  to  use  the  finite-element  method. 

The  recent  establishment  of  an  implicit  finite-element  solution  algorithm 
(Baker  and  Soliman,  1979,  1983)  has  substantially  advanced  the  solution 
technique.  Using  von  Neumann  linearized  stability  analysis,  the 
Baker-Soliman  algorithm  is  fourth-order  phase  accurate  with  third-order 
dissipation  in  the  large  Reynolds  number  range.  The  principle  of  the 
Baker-Soliman  (1979)  implicit  scheme  is  based  on  the  Galerkin  weighted 
residual  formulation.  For  parabolic  and  hyperbolic  partial  differential 
equations,  the  method  first  generates  linear,  quadratic,  and  two  cubic 
polynomials.  A  finite-difference  procedure  is  then  used  to  solve  the  resultant 
ordinary  differential  equation  system. 

From  a  theoretical  point  of  view,  modelers  who  use  the  finite-element 
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method,  emphasize  that  the  method  will  provide  a  thoroughly  structured 
procedure  to  transform  the  Navier  Stokes  equation  into  a  large  orderly 
equation  system  written  on  the  selected  discrete  variables  via  calculus  and 
vector-field  theory.  For  nonlinear  aerodynamic  problems,  the  finite  element 
formulation  is  usually  associated  with  the  inversion  of  large  matrices,  which  is 
deemed  to  be  computationally  undesirable.  Not  so  well  developed  in  the 
hypersonic  fluid-dynamic  modeling  field  as  the  finite-difference  method,  the 
finite  element  method  does  appear  to  give  robust  and  accurate  algorithms. 
Software  support  for  the  finite-element  models  are  relatively  limited  at  the 
present  time.  Recently  Bivens  (1989)  made  some  assessments  on  the 
reliability  of  using  finite  element  techniques. 


SPECTRAL  AND  PSEUDO-SPECTRAL  METHODS 

Relatively  speaking,  numerical  methods  such  as  the  finite-difference 
methods  (FDM)  or  the  finite-element  method  (FEM)  generally  use 
lower-order  difference  approximation  for  derivatives.  As  a  consequence,  phase 
errors  result.  Representing  the  derivatives  with  spectral  method  drastically 
reduces  these  kind  of  errors. 

Historically,  it  has  been  a  standard  technique  to  seek  a  solution  to  a 
dynamic  equation  as  a  series  of  known  functions,  so  that  each  spectral 
component  in  the  decomposition  of  the  solution  is  easier  to  solve  than  the 
complete  solution,  usually  via  an  analytical  method.  Laplace  and  Fourier 
transforms  were  often  used  to  reduce  the  degree  of  transcendency  for  solving 
numerous  differential  equations.  (Through  transformation,  partial  differential 
equations  become  ordinary  differential  equations,  for  example.) 

The  spectral  method  was  first  used  on  a  sphere  by  expansion  of  a  flow  field 
in  a  series  of  surface  harmonics  (Silberman,  1954).  The  main  difficulty  in 
employing  the  spectral  technique  at  that  time  seemed  to  be  the  considerable 
amount  of  arithmetic  operation,  and  the  computer  memory  required  to  store 
the  iteration  coefficients  for  each  mode.  In  the  pseudospectral  method 
(Orszag,  1970a),  the  mode  equations  are  not  actually  transformed  as  they  are 
in  the  spectral  method  (Orszag,  1970b),  but  the  derivatives  are  evaluated  by 
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the  fast  Fourier  transforms  (FFT).  The  FFT  method,  which  requires  only  N 
log^N  operations  per  time  step,  has  improved  the  operation  substantially  as 

compared  to  the  conventional  Fourier  transformation,  as  long  as  the  N  is  a 
power  of  2.  This  is  called  the  Cooley-Tukey  algorithm.  In  the  computation,  N 
is  the  number  of  modes  or  the  cut-off  frequency  to  be  considered. 
Theoretically,  spectral  simulation  is  infinite-order  more  accurate  than  the 
other  numerical  methods  in  the  aerodynamic  computation  because  errors 
decrease  more  rapidly  than  any  finite  power  of  1/N  as  the  N  goes  to  infinity. 
The  higher  the  value  of  N,  the  more  resolution  the  model  will  have,  but  at  the 
cost  of  higher  computational  effort,  however.  The  original  FFT  algorithm  of 
Cooley  and  Tukey  can  handle  series  only  with  the  cut-off  frequency  N  being 
the  power  of  2,  which  is  quite  inconvenient  in  some  cases.  More  recently,  with 
an  arbitrary-radix  algorithm  available,  the  amount  of  arithmetic  operation  is 
still  much  more  economical  than  the  conventional  Fourier  transform. 

For  aerodynamic  problems  with  periodic  boundary  conditions,  a  Fourier 
spectral  method  is  particularly  efficient  and  accurate.  It  can  also  be  used  in 
conjunction  with  multi-dimensional  Navier  Stokes  equations  using  spectral 
approximation  in  one  of  the  spatial  directions.  For  certain  applications,  this 
approach  has  been  shown  to  be  more  accurate  than  other  numerical  methods 
as  indicated  from  theoretical  analysis.  For  example,  when  computing  an 
aerodynamic  flow  field  passing  a  cylinder  or  hemispheric  cylinder,  Fourier 
spectral  approximations  can  be  used  to  compute  appropriate  derivatives  in  its 
circumferential  direction  and  a  finite-difference  scheme  in  the  other  directions. 
Using  a  spectral  approximation  of  the  derivatives  in  the  convective  terms,  it  is 
possible  to  achieve  a  particular  level  of  accuracy  with  fewer  circumferential 
points  than  is  possible  with  a  fourth-order  finite  difference  scheme  (Reddy, 
1983).  In  the  circumferential  direction,  the  components  of  derivatives  are 
approximated  spectrally  by  the  even  and  odd  discrete  Fourier  transforms. 

Let  f.  = 

where  0-  =  —  ,  j=0, 1,2,3, . N 

J  N 


-35- 


The  derivative  of  f  at  the  discrete  nodes  6-  can  be  approximated  by 

•J 


(M  =  T  n  a  cos  .  j  =  0,1,2 . N 

at'  J  n=l  ^  N 


where  the  spectral  coefficient  a^^  is 


'n  -  ^  j  I/f  °  ]  ’  “  -  .  N-1 


This  is  quite  similar  to  the  traditional  Fourier  transform  of  the  odd  (sin) 
function.  Sometimes  high-order  spectral  schemes  will  develop  nonlinear 
instabilities  in  strong  convective  air  flows.  Under  this  condition,  a  low-pass 
filter  has  to  be  used.  Gottlieb  et  al.  (1981)  proposed  a  filter  function  that 
operates  on  the  Fourier  coefficient  a^. 

The  spectral  method  can  also  be  used  in  conjunction  with  the  transformed 

computational  space  (i.e.,  For  example,  a  three-dimensional 

aerodynamic  flow  along  a  cylindrical  body  can  be  computed  using  a 

curvilinear  grid,  in  which  the  coordinate  ^  varies  streamwise  (along  the  body) 

and  T]  in  the  circumferential  direction.  When  a  derivative  of  a  component 
A  A 

estimator  (F)  is  calculated  in  the  transformed  direction  rj,  (ie,  dFldt]),  by  the 

spectral  method,  the  discrete  Fourier  transform  consists  of  the  usual  even  and 

odd  parts.  For  example,  in  modeling  the  aerodynamic  flow  field  over  a 

missile-shaped,  semi-sphere  cylinder,  Reddy  (1983)  used  even  functions  to 

represent  the  velocity  components  u,w  and  coordinates  x,z  in  the 

circumferential  direction  (i.e.,  ?/);  odd  functions  represent  v  and  y.  In  the 

stream— wise  direction  an  implicit  finite-difference  scheme  (Pulliam  and 

Steger,  1980)  is  used.  On  a  vector  computer,  the  spectral  method  requires 

approximately  15  percent  more  time  but  has  twice  the  accuracy  when  both 

were  compared  with  wind-tunnel  experiments  (Hsieh,  1976).  Therefore,  the 

overall  performance  of  this  spectral  method  was  concluded  to  be  40  percent 

more  efficient  than  a  fourth-order  finite-difference  scheme.  However,  it  is 
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this  writer’s  opinion  that  this  method  should  be  used  only  in  one  of  the 
curvilinear  coordinates;  to  quantify  the  overall  efficiency  may  not  be  so 
straightforward.  For  most  hypersomc  aerodynamic  applications,  the  spectral 
scheme  has  been  associated  with  problems  with  periodic  boundary  conditions. 

For  incompressible  subsonic  aerodynamic  boundary  layer  flow,  the 
pseudospectral  method  has  been  coupled  to  the  polynomial  subtraction 
technique  for  the  evaluation  of  spatial  derivatives  at  certain  nonperiodic 
boundary  problems  (Lee  and  Shi,  1986).  Their  technique  is  essentially  to 
subtract  the  polynomial  portion  of  the  solution  to  form  a  periodic  function. 
Potential  NASP  application  includes  the  circular-type  duct  flow  associated 
with  the  turbo-ramjet.  It  may  be  simulated  with  the  pseudospectral  method 
for  higher-order  accuracy  than  is  possible  with  the  traditional  finite-difference 
scheme,  which  is  usually  of  second  or  fourth  order  accuracy. 

In  terms  of  numerical  accuracy,  some  theoretical  comparisons  are  made 
between  the  spectral  method  and  the  finite-difference  methods.  For  the  sake 
of  convenience,  these  comparisons  are  usually  made  with  fluid  dynamic 
problems  in  which  accurate  analytical  solutions  exist.  One  such  problem  is 
the  aerodynamic  flow  around  a  sphere.  A  simple  case  is  that  the  flow  is  first 
subsonic,  then  supersonic,  and  finally  subsonic  again.  This  kind  of  flow  was 
first  analyzed  by  Ringleb  (Ringleb,  1940),  and  it  is  also  a  good  example  of  the 
shock-capturing  technique  used  in  the  spectral  method.  This  type  is 
presented  at  the  end  of  this  section  together  with  the  shock-capturing 
methods  associated  with  other  numerical  schemes.  Recently,  Canuto  et  al. 
(1988)  compared  solution  accuracy  between  the  spectral  method  and  the 
popular  second-order  MacCormack  finite  difference  scheme  for  the  simulation 
of  Ringleb  flow.  The  maximum  error  in  the  computed  pressure  against  the 
analytic  solution  for  the  transonic  and  supersonic  ranges  appears  in  Table  3. 
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Table  3 

Maximum  Error  in  Wall  Pressure  as  Computed  by 
Spectral  and  Second-Order  Finite-Difference  Solution 
(Adapted  from  Canuto  et  al.,  1988) 


Speed  range 

Grid 

MacCormack 

Spectral 

transonic 

8x4 

2.6  X  10~^ 

2.2  X  10"^ 

transonic 

16  X  8 

1.1  X  10"^ 

1.9  X  10"^ 

transonic 

32  X  16 

3.2  X  10"^ 

5.0  X  10"^ 

supersonic 

4x4 

2.2  X  10”^ 

7.5  X  10“^ 

supersonic 

8x8 

4.1  X  10"^ 

1.1  X  10"® 

supersonic 

16x  16 

1.0  X  10"^ 

6.6  X  10"^^ 

Ringleb  flow  was  selected  for  the  evaluation  of  numerical  accuracy  because 
an  exact  analytical  solution  exists  for  the  steady,  isentropic  flow  as  was 
explained  above. 

Even  though  spectral  methods  yield  more  information  about  the  exact 
solution  than  other  low-order  numerical  schemes,  when  the  exact  solution  is 
discontinuous  or  contains  large  gradients,  the  extra  information  is  then  hidden 
in  the  form  of  numerical  oscillations.  Typically,  when  the  spectral  method  is 
used  to  simulate  a  flow  with  shocks,  it  yields  an  oscillatory  solution.  The 
oscillations  occur  not  only  near  the  shock  but  all  over  the  field  because  the 
spectral  method  is  global  in  nature.  Like  artificial  dissipation  terms  in  the 
finite-difference  schemes,  diffusion  and  antidiffusion  terms  have  been  added  in 
the  spectral  methods  (Taylor  et  al.,  1981)  for  calculating  shock  waves, 
rarefaction  waves,  and  contact  surfaces.  This  point  is  discussed  further  below. 
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NUMERICAL  TREATMENT  OF  SHOCKS 

One  dynamic  characteristic  of  supersonic  flow  is  the  existence  of  shock 
waves.  To  resolve  shock  within  a  finite  computational  grid  structure  requires 
special  treatment.  These  methods  vary  according  to  the  basic  numerical 
solution  scheme  associated  with  the  simulation  model.  Several  techniques 
discussed  are  possible. 


Shock  Tracking  Technique 

In  the  traditional  numerical  analysis,  shock  surface  is  handled  like  an 
interior  boundary  of  discontinuity  (Courant  and  Friedrichs,  1948).  This 
method  is  more  conservative  in  terms  of  physical  state  but  it  also  has  some 
drawbacks  (Rizzi  and  Engquist,  1987).  One  disadvantage  is  that  the  nature 
and  interaction  of  discontinuity  have  to  be  known  ahead  of  time  to  set  up 
computational  points  in  the  model  to  handle  them  properly.  Because  of  the 
difficulties  in  programming,  the  method  has  been  replaced  by  the  other 
methods. 


The  Random  Choice  Method  of  Glimm 

In  solving  an  initial-value  problem  for  a  hyperbolic  system,  Glimm  (1965) 
used  a  constructive  random  choice  method  similar  to  one  developed  by 
Godunov  in  1959.  The  principle  is  to  make  a  change  in  the  grid  function 
averaged  for  a  certain  number  of  time  steps.  The  number  varies  according  to 
the  time  the  shock  passes  a  grid  point.  The  change  can  be  done  on  the 
average  time  steps  using  an  algorithm  based  on  random  numbers.  The 
method  had  worked  mainly  in  one-dimensional  problems  including  reacting 
gas  flow  (Chorin,  1977)  but  it  has  never  worked  well  for  multi-dimensional 
systems.  Chorin  also  worked  on  the  numerical  solution  of  Boltzmann’s 
equation  applied  to  the  problem  of  shock  structure  in  a  one-dimensional  flow. 
That  specific  numerical  method  is  more  accurate  than  the  Monte  Carlo 
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methods  for  solving  Boltzmaim’s  equation  particularly  in  the  high  Mach 
number  region  (Chorin,  1972). 


Shock-Capturing  Technique 


This  method  is  used  in  the  majority  of  numerical  models  in  the  hypersonic 
aerodynamic  simulation  associated  with  NASP  analysis.  Its  major  advantage 
is  its  simplicity  in  programming,  which  assumes  no  prior  knowledge  of  the 
nature  and  the  location  of  the  shock  front.  The  shock-capturing  technique 
was  originally  proposed  by  von  Neumann  during  the  1940s.  Substantial 
progress  has  been  made  during  the  last  two  decades  (Godunov,  1970;  Jameson, 
1989).  The  method  is  based  on  the  principle  that  we  have  to  make  a 
compromise  between  the  accuracy  of  the  solution  value  and  the  accuracy  of 
the  shock  location.  The  final  result  is  to  give  up  a  certain  degree  (first  order) 
of  accuracy  in  the  immediate  vicinity  of  the  sharp  discontinuity.  However,  if 
the  numerical  scheme  used  for  simulating  hypersonic  flow  is  nondissipative, 
numerical  \iscosity  must  be  artificially  added  to  capture  the  shock  wave. 

To  add  artificial  viscosity,  it  is  desirable  that  the  shock  transition  layer 
not  extend  over  more  than  a  few  grid  spadngs  and  that  the  artificial  viscosity 
be  independent  of  the  shock  strength.  It  is  also  desirable  that  the  transition 
layer  travel  at  very  nearly  the  correct  speed  through  the  air.  Sometimes  it  is 
also  desirable  to  perform  some  numerical  filtering  to  smooth  out  short  waves 
caused  by  nonlinear  instabilities  in  the  smooth  air  flow  away  from  the  shock. 
The  numerical  filter  will  be  turned  off  automaticcdly  when  the  shock  wave  is 
not  present.  In  a  numerical  study  of  shock  boundary  layer  interaction, 
MacCormack  and  Baldwin  (1975)  devised  a  normalized  second-difference 
sensor  using  the  computed  pressure  across  three  grid  spaces  of  the  form: 


‘  I  Pj+1  +  2Pj  +  Pj.ll 


(  3.18  ) 


For  applications  in  two-dimensional  and  three-dimensional  supersonic 
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flows,  the  same  principle  can  also  be  extended  and  is  carried  out  dimension  by 
dimension.  The  same  concept  can  be  implemented  in  models  using  total 
variation  diminishing  (TDV)  schemes.  General  conditions  on  coefficients  that 
result  in  total  variation  diminishing  were  proved  by  Jameson  and  Lax  (1984) 
(see  also  Jameson,  1989;  Yee,  1987,  1989a  and  1989b;  Yee  and  Shinn,  1986; 
Yee  et  al.,  1988,  1990).  TDV  methods  can  give  sharp  discrete  shock  waves 
without  oscillations,  but,  as  expected,  the  price  for  obtaining  a  stable  solution 
by  filtering  is  the  slight  loss  of  accuracy  of  the  numerical  solution. 

Capturing  vortex  sheets  in  a  numerical  simulation  is  somewhat  more 
involved.  Presently,  no  well-estabUshed  technique  seems  to  be  available. 


Spatial  Switching  for  Shock  Resolution  in  the  ADI  Scheme 

The  class  of  ADI  schemes  discussed  in  the  previous  subsection  is 
nondissipative,  so  that  they  are  not  good  for  hyperbolic  equations  when  shock 
waves  occur  as  in  hypersonic  aerod)mamic  conditions.  Some  dissipative  terms 
have  to  be  added.  One  of  the  often-ased  algorithms  is  to  let 
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According  to  the  von  Neumann  stability  condition,  a  one-dimensional 
implicit  numerical  scheme  is  stable  if  the  value  of  ui  is  between  zero  and  unity. 
For  multi-dimensional  cases,  some  experiments  have  to  be  made.  The 
often-used  value  of  u  has  been  0.5. 

Adding  a  dissipative  term  to  reduce  post-shock  oscillations  works  under 
certain  situations.  But  there  are  cases  in  which  adding  a  dissipative  term  has 
to  be  accompanied  by  another  treatment  for  resolving  shock  waves.  Sometime 
a  second-order  explicit  step  which  switches  difference  operator  from  "central" 
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to  "upwind"  (one-sided)  across  a  discontinuity  can  greatly  reduce  the  spurious 
oscillations.  The  method  is  quite  commonly  adapted  in  CFD  modeling  of 
marine  pollutant  transport  simulation  if  a  pollutant  with  high  concentration  is 
injected  into  currents. 


Shock-Capturing  and  Fitting  Technique  in  Spectral  Method 

Even  though  the  spectral  method  generally  gives  more  accurate  numerical 
results,  the  ability  for  resolving  the  shock  is  generally  not  as  good  as  the 
time-domain  numerical  schemes  such  as  the  finite-difference  method.  In 
spectral  simulations  shock  waves  are  usually  represented  as  a  series  of 
osciUations.  To  resolve  the  shock,  different  types  of  smoothing  filters  are 
usually  employed  as  a  post-processing  step.  Commonly  used  filters  include 
exponential,  Lanczo,  or  other  types  of  cosine  filters  (Canute  et  al.,  1988). 
Numerical  filtering  across  the  frequency  domain  between  0  and  2ir  usually 
resolves  the  location  of  the  shock. 

Another  technique  to  solve  the  shock-induced  oscillation  in  the  spectral 
method  is  to  treat  the  shock  front  as  an  internal  computational  boundary. 

The  shape  of  the  shock  front  is  determined  during  the  computation.  This 
technique  is  termed  shock-fitting  (Moretti,  1968,  1972;  Salas  et  al.,  1982; 

Zang  et  al.,  1984;  Hussaini  et  al.,  1985;  Canute  et  al.,  1988).  Up  to  now,  the 
spectral  shock-fitting  schemes  have  been  applied  only  to  the  numerical  models 
based  on  the  Euler  equations,  which  neglect  viscous  terms. 


Numerical  Treatment  of  Subsonic  Regions  in  a  Supersonic  Flow 

When  an  aerospace  plane  flies  at  supersonic  speed,  its  configuration  may 
consist  of  some  subsonic  regions  (pockets)  in  the  generally  supersonic  flow. 
These  areas  could  be  developed  as  a  result  of  gradual  compression  or  caused 
by  flow  separation.  These  phenomena  need  special  numerical  treatments. 
The  traditional  approach  is  to  use  different  marching  codes  for  space  and  for 
time.  Recently,  Chakravarthy  et  al.  (1988)  handled  this  problem  by  using  a 
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unified  marching  scheme  for  both  space  and  time  using  the  relaxation  method. 
In  their  unified  solution  algorithm  (USA  code),  the  FNS  solver  is  based  on  the 
total  variation  diminishing  (TVD)  formulation  using  the  finite-volume 
approach.  With  the  Baldwin/Lomax  and  one-equation  turbulence  closure 
method  (see  Sec.5),  a  fully  upwind,  not-flux-limited  scheme  is  used  in  the 
supersonic  region  external  to  the  subsonic  part  of  the  boundary  layer.  This 
region  involves  only  a  forward  marching  sweep.  For  subsonic  regions,  the 
solution  is  to  first  march  forward  with  one  or  two  subiterations,  then  to  follow 
with  a  backward  marching  sweep.  The  TVD  formulation  as  described  above 
captures  shock  wave  in  the  subsonic  region. 
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4.  GRID  GENERATION  AND  COORDINATE  TRANSFORMATION 

TECHNIQUES 


INTRODUCTION 

One  major  advance  in  the  field  of  aerodynamic  modeling  has  been  the 
development  of  a  grid-generation  technique.  This  technique,  when  combined 
with  efficient  finite-difference  methods,  gradually  evolves  into  a  major 
aircraft  design  tool  supplementary  to  ground  experiments.  The  difference 
between  a  grid-generation  system  and  a  moving  adaptive  grid-generation 
system  is  that  the  former  generates  an  aerodynamic  computational  network 
before  carrying  out  the  calculation,  the  latter  generates  grid  during 
computation  according  to  certain  criteria  such  as  the  pressure  gradient. 
However,  the  adaptive  methods  often  use  the  regular  grid-generation  method 
to  create  the  initial  grid  as  the  starting  condition.  At  present,  nearly  all  the 
NASP-related  hypersonic  aerodynamic  computations  involve  grid  generation 
of  some  type.  In  many  cases,  the  grid  system  involves  the  segmentation  of 
the  aerodynamic  flow  field  into  subregions,  with  grids  being  generated  in  each 
subregion.  Continuity  is  then  enforced  at  the  interface.  This  section  reviews 
the  methods  commonly  used  in  numerical  grid  generation  for  the  NASP 
applications. 


METHODS  AND  PRINCIPLES  OF  GRID  GENERATION 

An  accurate  description  of  the  numerical  grid  generation  process, 
according  to  Thompson  (1984),  may  be  described  as  follows: 

Difference  representations  on  curvilinear  coordinate  systems  are 
constructed  by  first  transforming  derivatives  with  respext  to 
Cartesian  coordinates  into  expressions  involving  derivatives  with 
respect  to  the  curvilinear  coordinates  and  derivatives  of  the  Cartesian 
coordinates  with  respect  to  the  curvilinear  (metric)  coefficients.  The 
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derivatives  with  respect  to  the  curvilinear  coordinates  are  then  replaced 
with  difference  expressions  on  the  uniform  grid  in  the  transformed 
region. 

A  grid  system  can  usually  be  generated  using  one  of  the  following 
methods  (Thompson,  1984;  Thompson  et  al.,  1985,  Thompson  and  Ferziger, 
1989); 


•  Algebraic  system, 

•  Conformal  mapping  system, 

•  Elliptic,  parabolic  and  hyperbolic  systems. 

Three-dimensional,  mass-averaged  Navier  Stokes  equations  in  a 
transformed  space  are  presented  in  Sec.  2. 


ALGEBRAIC  SYSTEM 

In  an  algebraically  generated  grid  system,  the  spacings  are  basically 
interpolated  among  boundaries.  This  type  of  generating  method  is 
well-suited  when  used  together  with  computer-aided-design 
(CAD/interactive  graphic  system).  The  characteristics  of  this  method  are: 


•  Fast, 

•  Allows  explicit  control  of  grid  distribution,  or 

•  No  inherent  smoothing  mechanism  (even  though  cubic  spline  may  be 

used), 

•  Often  uses  transfinite  interpolation  as  for  sculptured  surfaces.  This 

is  similar  as  in  the  CAD  system  which  primarily  uses  the  Boolean  sum 
projector  method  from  the  boundaries. 

In  one-dimensional  cases,  the  interpolation  scheme  often  included 
Lagrangean  polynomial  interpolation,  Hermitian  interpolation,  cubic  spline, 
tension  spline,  B-spline,  etc.  For  two-  to  three-dimensional  cases,  transfinite 
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interpolations  are  usually  applied,  for  example,  if  linear  interpolation 
functions  are  needed  in  each  curvilinear  direction  (2-D,  N  by  M)  from  a  set  N 
+  M  intersecting  curves.  These  functions  are: 


(4.1) 


(4.2) 


The  interpolation  matches  the  function  on  the  boundary  defined  by  ^=0 
and  ^=I  in  the  first  equation,  or  j/=0  and  7?=J  in  the  second  equation.  The 
interpolation  is  transfinite  if  the  match  on  the  entire  boundary  at  a 
non-denumerable  number  of  points. 

The  general  form  of  the  transfinite  interpolation,  which  gives  algebraically 
the  best  approximation,  is; 

N  M 

=  I  (-f-)  +  I  l'’m  (y) 

n=l  m=l 

N  M 

-  I  I  (^-3) 

n=lm=l 

or 


(4.4) 


The  first  term  is  the  result  at  each  point  in  the  field  of  the  unidirectional 
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interpolation  in  the  7/-(iirection  and  the  bracket  is  the  difference  between  the 
specified  values  on  the  lines  and  the  result  of  the  unidirectional 

interpolation  on  these  lines.  The  spline-blended  norm  gives  the  smoothest 
grid  with  continuous  second  derivatives  according  to  Thompson  (1984). 


CONFORMAL  MAPPING  TECHNIQUE 

Conformal  mapping  is  a  classic  technique  based  on  the  principle  of 
transformation  between  two  regions  that  are  conformally  (one-to-one) 
equivalent.  This  class  of  methods  includes  the  solution  of  integral  equations, 
the  expansion  of  power  series  or  Fourier  series,  and  the  construction  of 
Schwarz-Christoffel  transformation  (see,  e.g.,  Sokolinkoff  and  Redheffer, 
1966).  Conformal  systems  have  the  advantage  of  introducing  the  fewest 
additional  terms  in  the  transformed  partial  differential  equations.  Even 
though  systems  generated  by  conformal  mapping  are  inherently 
two-dimensional,  more  complicated  shapes  have  recently  been  constructed 
(Thompson  and  Ferziger,  1989).  Conformal  mappings  do  not  exist  in  three 
dimensions  except  in  trivial  cases.  A  curvilinear  coordinate  system  generated 
by  a  conformal  mapping  is  usually  quite  rigid  in  a  way  that  little  control  can 
be  exerted  over  the  distribution  of  grid  network. 

Aerodynamic  shapes  can  sometimes  be  transformed  into  analytical  forms 
such  as  a  circle.  For  example,  an  airfoil  shape  can  be  approximated  as  the 
image  of  a  circle  through  the  transformation: 

z  =  ^  +  -^  (4.5) 

Under  the  inverse  transformation,  a  given  airfoil  will  map  to  a  curve  that 
is  nearly  circular.  Analytical  functions  are  also  useful  to  generate  grid  near 
the  boundaries  with  slope  discontinuities.  If  algebraic  methods  are  used,  the 
discontinuities  will  propagate  into  the  physical  region,  thus  inducing 
nonsmooth  grid  spacings. 

Complex  aerodynamic  shapes  can  be  treated  by  a  sequence  of 


transformations  each  to  an  analytical  function  in  succession  or  in  isolation 
through  an  iterative  procedure  (Thompson  et  al.,  1985). 

ELLIPTIC,  PARABOLIC,  AND  HYPERBOLIC  SYSTEMS 

These  systems  are  classified  according  to  the  characteristics  of  the 
generating  equation.  For  example,  Laplace  and  Poisson  partial  differential 
equations  are  of  the  elliptic  type.  This  type  of  system  is  the  most  frquently 
applied  among  the  three.  It  has  many  convenient  features.  One  of  the  most 
important  characteristics  of  the  Laplace  systems  is  that  it  guarantees  a 
one-to-one  mapping  for  boundary-conforming  curvilinear  coordinate  systems 
on  generally  closed  boundaries.  It  usually  generates  a  very  smooth  grid 
network.  The  difference  between  a  Laplace  and  a  Poisson  system  is  the 


generating  equation: 

Laplace  equation 

O 

II 

0=1,2, 3) 

(4.6) 

Poisson  equation 

=  P* 

(i=l,2,3) 

(4.7) 

P^  in  Eq.  4.7  are  called  the  control  functions,  which  are  used  to  control 
the  spacing  and  the  direction  of  the  coordinate  lines.  Since  the  Laplacian 
operator  lacks  control  functions,  the  coordinate  lines  will  tend  to  be  generally 
equally  spaced  away  from  the  solid  model  boundaries  near  the  far  field. 
Sometime  this  is  a  desirable  property.  The  selection  of  control  functions  for 
the  Poisson  system  is  described  in  Thompson  et  al.  (1985). 

Grid-generating  systems  can  also  be  based  on  the  numerical  solutions  of 
parabolic  and  hyperbolic  partial  differential  equations.  In  the  solution 
sequence,  one  proceeds  in  the  direction  of  one  curvilinear  coordinate  between 
two  boundaries  for  the  two-dimensional  case  and  between  the  two  boundary 
surfaces  for  a  three-dimensional  case. 

In  elliptic  grid  generating  systems  as  shown  above,  second  derivatives 
exist  in  both  directions  (as  indicated  by  the  operator).  Therefore,  only  in 
an  elliptic  system  can  the  entire  boundaries  of  a  general  region  be  specified. 


The  hyperbolic  system,  on  the  other  hand,  allows  only  one  boundary  to  be 
specified.  And  because  of  this,  the  hyperbolic  system  is  faster  than  the 
elliptic  system  by  one  or  two  orders  of  magnitude.  The  same  is  true  for  a 
parabolic  system  in  which  a  second  derivative  in  one  of  the  directions  does  not 
appear  (like  the  heat  equation).  Parabolic  and  hyperbolic  systems  are  similar. 
One  major  difference  is  that  in  a  parabolic  system,  influence  from  the  other 
boundary  still  exists  in  the  governing  equation.  The  speed  advantage  of  the 
parabolic  system  over  the  elliptic  system  is  roughly  the  same  as  that  of  the 
hyperbolic  system.  Marching  grid  generation  using  hyperbolic  partial 
differential  equations  was  described  in  some  detail  by  Steger  et  al.  (1977). 

The  Poisson  generating  system  is  perhaps  the  most  applied  elliptic  system. 
We  will  discuss  it  here  in  some  detail.  In  a  Poisson  (inhomogeneous  elliptic) 
system  (Eq.  4.7) 


=  pi  (i= 1,2,3) 

The  control  functions  P*  are  used  to  control  the  spacing  and  direction  of 
the  coordinate  lines.  The  selection  of  a  control  function  for  the  Poisson 
system  where  P*  is 


jSl  kSl 


(4.8) 


where  g^  are  the  components  of  the  contravariant  metric  tensor  which  are 
the  dot  products  of  the  contravariant  base  vectors  (a')  of  the  curvilinear 

coordinate  system. 


(4.9) 


Based  on  this  relationship,  the  Poisson  generation  system  can  be  defined 
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3  3 


=  I  I  ^ik  (  ^ 

j4l  k=l  ■ 


(4.10) 


When  the  computation  is  carried  out  in  a  rectangular  transform  space,  the 
curvilinear  coordinates  are  the  independent  variables  with  the  cartesian 
coordinates  x-  as  dependent  variables.  The  more  common  form  used  as  the 

generating  system  in  the  transformed  region  (Thompson  et  al.,  1985)  with 
three  control  functions  is 


3  3  3 

I  I  1  e%'k  =  ‘> 

iAi  jsi  -{"{J  k=i  *  'r 


(4.11) 


The  two-dimensional  form  of  the  generation  system  with  two  control 
functions  is 


+  P  ??)  +  8u(!„+  Q  £,)  -  If,  =  « 


(4.12) 


In  the  physical  space,  the  system  is 


V^«  =  -^P 


(4.13) 


Q 

6 


(4.14) 


For  the  orthogonal  system  (612=613=123=0.  611=622=633=““®*  ). 

12 

two-dimensional  (g22=l),  in  physical  space:(g  =0) 


"lr(  h^  M )  + 


(4.15) 
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where 


h 

h 


dr\  \  ,  d  /  ^2 

-^>  =  0 

(4.16) 

1=  /g]7 

(4.17) 

2  “  ^  ^22 

(4.18) 

gjj  are  elements  of  the  metric  tensor.  The  ratio  ^l/h2  and  ^2/hj  are  the 
attraction  functions  that  need  to  be  prescribed.  These  ratios  become 

h,  _ ___ 

~E~  ^  811/822  =  ^ 

2 

A  curvilinear  grid  network  can  be  generated  by  first  establishing  a  coarse 
grid  using  a  graphic  input  device  such  as  a  CAD  system.  The  ratio  of  A^/At; 
over  the  computational  field  can  be  interpolated  into  a  finer  ratio  where  ^  and 
Tj  have  to  be  computed.  The  attraction  function  q  is  calculated  using  Eq. 

4.19. 

Substituting  a  and  1/a  for  hj/h2  and  h2/h2  in  Eqs.  (4.15)  and  (4.16),  an 

orthogonal  curvilinear  computational  grid  is  then  generated. 

Since  the  Jacobian  of  the  transformation  for  a  two-dimensional 
orthogonal  system  is 


J  =  /  gi"  =  V  giig22  =  ^ 

2i0.5 


=  (- 


dx  dy  ' 

r  dx  ] 

hr\ 

+ 


(4.21) 


The  continuity  equation  for  a  two-^iimensional  system  with  free  surface  is 
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-I-  + +;|^-|-[(<l+0v/6^]  =0  (4.22) 

where  the  transformation  coefficients  and  usuaUy  located  at  the 

V  and  u  vectors,  respectively,  to  maintain  the  conservation  of  mass  during  the 
integration  process.  Within  a  control  volume,  they  are  located  at 


QU,v,^,y) 


(4.23) 

(4.24) 


NUMERICAL  SOLUTION  SCHEME  -  POISSON  SYSTEM 

The  solution  of  an  elliptic  Poisson  equation  is  a  boundary  value  problem. 

If  ^  is  a  prescribed  function  on  the  boundary,  such  a  problem  is  called  a 

"Dirichlet  problem."  If,  instead  of  the  value  the  value  of  dijdv,  the  normal 
derivative  of  is  prescribed  on  the  boundary,  the  problem  is  called  a 
"Neumann  problem."  Several  numerical  iteration  methods  such  as  the 
Richardson  and  Liebmann  line  iteration  scheme  can  be  used.  The  method 
developed  by  Peaceman  and  Rachford  (1955)  seems  to  be  the  quickest 
iterative  method  in  which  line  iteration  schemes  are  used  in  the  columns  and 
rows  alternatively.  The  explicit  description  of  the  method  is  contained  in  the 
following  recursion  formula: 
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=  -  ei,j-i''“>  +  (2-^)?i/“)-  {i,j+/"”>+  ^  fy  (4-25) 

=  -  +  (2-/.n)«i_/^"+l>-  ei+ij<^"+'>+  ^  fy  (4.26) 


where  n=0,l,2,...,  and  hopefully  the  sequence  converges.  The  value  ^  is  an 
extrapolation  parameter  that  is  to  be  determined  so  that  the  method  will 
converge  as  quickly  as  possible.  Peaceman  and  Rachford  suggest  putting  pn 
=pk  if  n  E  k  (mod  p)  where 


pk  =  4  sin^-^^^^p 


(4.27) 


ORTHOGONAL  SYSTEMS 

Orthogonal  systems  are  important  to  grid  generation  for  two  major 
reasons.  First,  the  truncation  error  associated  with  an  orthogonal  system  in 
the  difference  expression  is  minimal.  Therefore,  the  method  is  more  accurate 
than  others.  Second,  the  orthogonal  system  induces  fewer  additional  terms  in 
the  transformed  partial  equations  to  account  for  the  effects  of  curvature  and 
centrifugal  forces.  Therefore,  the  method  is  computationally  very  efficient. 

As  a  consequence,  perhaps  more  aerodynamic  applications  are  based  on 
orthogonality  than  on  any  other  method. 

A  true  orthogonal  system  in  a  three-dimensional  case  is  also  very  difficult 
to  implement.  For  aerospace  applications,  orthogonality  will  likely  be  over  a 
surface  coordinate  network. 
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MOVING  ADAPTIVE  GRID  SYSTEMS 

One  of  the  newest  developments  in  the  field  of  hypersonic  aerodynamic 
modeling  is  in  the  technique  of  moving  adaptive  grid  generation.  During 
simulation,  the  grid  network  moves  as  the  physical  solution  develops.  By 
sensing  the  gradient  of  the  solution,  the  grid  concentrates  in  regions  of  large 
variations. 

The  grid  adaptation  can  generaUy  be  achieved  by  arranging  the 
neighboring  points  in  an  identifiable  order  so  that  a  continuous  function  can 
be  represented,  and  errors  evaluated  and  redistributed  while  the  grid  network 
is  shifting.  In  doing  this,  however,  sufficient  resolution  in  time  and  space  is 
required  to  minimize  and  to  evaluate  truncation  error  so  that  numerical 
oscillation  can  be  avoided. 

In  an  economical  sense,  the  method  offers  better  resolution  at  fewer  points 
than  the  other  nonadaptive  schemes,  even  though  somewhat  more 
computations  are  involved.  With  the  potential  saving  in  the  total  number  of 
points  for  a  given  level  of  accuracy,  the  adaptive  method  is  still  competitive 
with  the  other  schemes  in  terms  of  computer  time.  The  method  is  most 
suitable  for  problems  for  which  we  have  no  prior  knowledge  of  their  solution. 
However,  an  initial  fixed  grid  has  to  be  generated  with  the  other  methods. 

In  solving  multi-dimensional  aerodynamic  problems  such  as  those 
associated  with  NASP,  if  the  variability  of  the  solution  sought  is  mainly  in 
one  direction,  then  the  grid  adaption  can  be  applied  with  the  grid  location 
constrained  to  move  along  one  direction  only.  The  grid  spacing  along  the 
region  of  large  gradients  is  placed  according  to  the  principle  of 
equi-distribution  of  certain  weighting  functions.  The  criterion  is  to  select 
these  weighting  functions  to  be  inversely  proportional  to  the  grid  spacing. 
Therefore,  larger  gradients  are  resolved  by  smaller  spacing.  We  will  use  this 
rather  simple  case  for  illustration  purpose. 

Let  x(i)  represent  grid  spacing  and  Wj  denote  a  weighting  function.  The 

idea  is  that  the  product  of  the  weighting  function  and  the  grid  spacing  would 
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be  a  constant  (Thompson  et  al.,  1985).  In  the  transformed  domain, 


A^.w 


•w  =  const. 


(4.28) 


where  is  grid  spacing 

If  a  velocity  gradient  is  selected  as  the  weighting  function,  it  can  be 
selected  as 


w  = 


da 

W'  J 


(4.29) 


The  above  formulation  has  the  advantage  that  the  resulting  spacing  near 
the  far-Celd  boundary  will  not  be  very  large.  Ideally,  it  would  be  more 
desirable  to  have  equal  spacing  if  the  gradient  is  zero  such  as  near  the  far 
field.  To  achieve  this,  one  has  to  select  a  set  of  parameters  according  to  the 
curvature  of  the  solution  gradient  (Eiseman,  1985)  such  as: 


w 


1  + 


^2|k|  [l-Ka2(5u/dx)2j 


1/2 


(4.30) 


where  a  and  0  are  parameters  to  be  selected  and  k  is  the  curvature  of  the 
solution  curve. 

Three-dimensional  adaptation  is  much  more  complicated  than  the  above 
example.  However,  the  general  idea  is  still  the  same.  Recently,  tension 
spring  analog  has  been  used  for  the  three-dimensional  adaptation  (Nakahashi 
and  Deiwert,  1986)  in  which  the  weighting  function  is  treated  as  the  tension 
spring  constant  of  the  three-dimensional  spring  system.  In  the  system,  each 
grid  node  is  suspended  by  six  tension  and  twelve  torsion  springs.  Rather  than 
solving  it  simultaneously,  the  system  splits  into  a  sequence  of 
one-dimensional  adaptations  in  which  three-dimensional  grid  movements  are 
achieved  by  successive  applications  of  the  one-dimensional  method.  The 
weighting  function  w  (i.e.,  the  tension  spring  constant)  is  derived  according  to 
the  solution  gradient  over  an  arc  length  multiplied  by  a  group  of  tension  and 
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torsion  spring  coefficients.  The  tension  spring  coefficient  affects  the  mesh 
spacing  along  each  line.  It  is  selected  by  considering  the  range  of  spacing 
change  (max-min)  along  a  line.  Torsion  spring  constants  control  the  damping 
between  tension  and  torsion  forces.  As  a  result,  larger  values  give  a  more 
gradual  change  in  the  mesh  spacing. 


ADAPTIVE  MOVING  FINITE  ELEMENTS 

This  dynamically  adaptive  method,  in  essence,  adds  dependent  variables 
at  the  location  where  the  Galerkin  weighted  residual  process  is  used.  In  the 
aerodynamic  flow  formulation,  the  entire  domain  is  reduced  into  finite 
elements  whose  contributing  residuals  are  evaluated  by  trial  solution.  At 
each  grid  point  and  on  each  element,  the  residual  is  required  to  be  orthogonal 
to  all  the  basic  functions.  As  a  result,  the  location  of  grid  points  becomes 
part  of  the  solution.  At  the  present  time,  solution  methods  are  developed 
only  for  two-dimensional  systems.  For  more  complex  systems,  it  seems 
apparent  that  roundoff  error  can  play  a  major  part  in  the  accuracy  of  the 
result,  inasmuch  as  the  complexity  of  a  large  problem  necessitates  a  larger 
number  of  iterations  both  within  and  per  time  step. 

For  NASP  application  with  shocks,  solution  derivatives  of  the  density  field 

seem  to  be  a  good  choice  for  the  weighting  function.  Gnoffo  (1983)  used 

Mach  numbers  for  this  purpose  to  compute  the  flow  field  over  the  Galileo 

5  5 

probe  where  values  were  computed  for  M=50,  Re=10  -5x10  (7=1.4, 

Pr=0.72,  Sotherland’s  law  for  viscosity).  Experiment  verification  for  M=6 
(Libby  and  Cresd,  1961)  was  made.  For  the  hypersonic  flow  simulation,  the 
adaptation  technique  is  based  also  on  an  equivalent  spring  analogy  where 
springs  connect  adjacent  mesh  points  and  spring  constants  are  a  function  of 
the  user-specified  gradient  between  the  point.  The  scheme  is  coupled  to  the 
finite  volume  method.  In  the  formulation,  the  integral  form  of  the  governing 
conservative  laws  is  approximated  on  cells  whose  corners  are  defined  by  the 
position  of  grid  points  in  physical  space.  On  a  rectangular  grid,  the  algorithm 
reduced  exactly  to  the  MacCormack’s  explicit  method  (Sec.  3)  which  is 
second-order  accurate  both  in  time  and  space.  During  adaptation,  certain 
filtering  of  the  spring  constants  was  involved. 
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Recently,  Peraire  et  al.  (1987)  introduced  an  adaptive  remeshing 
procedure  for  improving  the  quality  of  steady-state  solutions  to  the 
two-dimensional  Euler  equation  with  a  finite-element  algorithm.  The 
method  involves  an  iterative  step  using  the  computed  solution  to  determine 
the  optimal  value  according  to  an  indicator  of  the  error  magnitude  and 
direction.  The  scheme  has  the  advantage  of  gradually  improving  the  quality 
of  the  solution  without  significantly  increasing  the  total  number  of  unknowns 
at  each  stage  of  iteration.  No  results  are  available  from  the  three-dimensional 
analysis  yet. 
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5.  MODELING  TURBULENCE 


INTRODUCTION 

One  major  effort  in  hypersonic  aircraft  design  is  to  predict  the  location  of 
laminar-turbulent  transition.  The  range  and  the  magnitude  of  uncertainty  in 
the  prediction  increase  as  the  Mach  number  increases.  The  uncertainty 
translates  directly  into  uncertainties  in  the  prediction  of  vehicle  performance, 
weight,  and  other  control  parameters.  To  narrow  the  range  of  uncertainty  in 
transition  prediction,  the  numerical  modeling  of  turbulent  transitional 
processes  is  extremely  important. 

The  National  Research  Council  made  an  objective  assessment  of  the 
capabilities  and  future  directions  in  CFD  (NRC,  1986).  In  the  area  of 
modeling  turbulence,  they  have  reached  the  following  conclusion: 

Turbulence  modeling,  including  modeling  of  the  laminar-turbulent 
transition,  is  becoming  a  padng  technology.  Present  turbulent 
models  are  adquate  only  for  use  in  relatively  simple  flows,  and  do 
poorly  in  flows  with  strong  three-dimensionality,  massive  separation, 
large-scale  unsteadiness,  strong  density  gradients,  strong  rotation, 
and  chemical  reaction.  Large  eddy  simulations  have  not  yet  been 
developed  for  boundary  conditions,  geometries,  and  Mach  numbers  of 
practical  interest.  (NRC,  1986) 

We  concur  with  the  NRC  and  observe  that  additional  research  in  the 
modeling  of  turbulence  is  needed  to  reduce  the  uncertainties  in  predicting  the 
performance  parameters  of  a  vehicle  if  CFD  simulation  is  used  as  a  design 
tool.  In  this  section,  we  will  examine  various  turbulence  models  used  in 
hypersonic  research  and  try  to  make  some  quantitative  assessment  of 
uncertainties  associated  with  each  approach. 
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A  fflSTORICAL  BACKGROUND  OF  MODEL  TURBULENCE 

The  original  Navier  Stokes  equation  (Navier,  1822;  Poisson,  1829; 
Saint-Venant,  1843;  and  Stokes,  1845)  was  formulated  only  for  describing  the 
laminar  flow  field  in  which  the  momentum  transfer  due  to  viscous  effects  is 
included  in  the  viscous  terms.  The  amount  of  momentum  loss  or  gain  due  to 
the  viscous  effect  is  proportional  to  the  coefficient  of  viscosity  (/z)  of  the  fluid 
[gr  cm'^sec”^]  which  is  assumed  constant.  The  equations  are  (see  Sec.  2  for 
a  generalization): 


d[x  ,  ..  d\i  ,  ^  „  d\i 

-5r  +  "-ffir  + +  *-3r 


dw  ,  „  dv  dv  ,  dv 

~dr  + 


dw  dw  ,  dw  ,  dw 


dw 

w 


■^-^  +  -^V2v 

fi  oy  p 

■  J _ iE_+ JLc!* 

p  p  * 


(5.1) 

(5.2) 

(5.3) 


The  last  term  in  each  equation  contains  the  kinematic  (molecular) 
viscosity.  The  Navier  Stokes  equations,  including  mass,  energy,  and  species 
conservation,  can  also  be  derived  from  the  kinetic  theory  of  gases  using  the 
Chapman  and  Enskog  procedure  (Chapman  and  Cowling,  1939;  Tsien,  1958). 
Presumably,  they  apply  to  reacting  gas  flows,  both  turbulent  and  laminar, 
under  compressible  flow  conditions.  For  most  aerodynamic  modeling 
applications,  the  molecular  viscosity  of  single  component  gas  is  often 
calculated  using  Southerland’s  law  with  appropriate  constants  for  the  species 
under  consideration.  For  the  case  of  multiple  component  gases  such  as  the 
mixture  of  hydrogen  and  air,  its  binary  molecular  viscosity  can  be  calculated 
by  the  Wilke’s  formula. 

Turbulent  flows,  on  the  other  hand,  cannot  be  computed  with  these 
time-dependent  Navier  Stokes  equations  unless  the  flow  field  is  resolved  by  a 
numerical  model  whose  computational  network  approaches  Re®/^  grid  points 
(Case  et  al.,  1973).  With  that  many  points,  the  effects  of  even  the  smallest 
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turbulent  eddies  can  then  be  represented  entirely  by  these  viscosity  terms. 

For  turbulent  flow,  it  was  Boussinesq  (1877)  who  used  the  "eddy-viscosity" 
concept  which  assumes  that  the  gross  effects  of  both  the  turbulent  and  the 
laminar  portions  of  the  viscous  effects  are  proportional  to  the  mean  velocity 
gradient.  Eddy  viscosity  is  analogous  to  the  viscous  stress  coefficient  (/i) 

in  the  laminar  flow.  Although  /i  is  a  fluid  property,  eddy  viscosity  depends  on 
the  velocity  gradients  vnthin  the  flow.  However,  Boussinesq’s  eddy  viscosity 
concept  forms  the  foundation  of  turbulence  theory.  Together  with  Prandtl’s 
contribution,  the  eddy  viscosity  type  of  closure  methods  such  as  the 
zero-equation,  one-equation  and  two-equation  turbulence  models  are  the 
primary  method  used  in  the  aerodynamic  flow  modeling.  In  these  models,  the 
shear  stresses  in  the  Navier  Stokes  equations  are  computed  as  a  function  of 
mean  velocity  gradient  and  the  eddy  viscosity,  of  the  general  form: 

Since  the  eddy  viscosity  varies  within  the  flow  field  under  different 
situations,  these  zero-,  one-,  and  two-equation  models  are  designed  to 
account  for  the  variability  of  these  coefficients,  using  hypotheses, 
assumptions,  and  experimental  results  in  the  process  of  seeking  certain 
universal  relationships.  However,  the  ultimate  solution  to  account  for  the 
various  flow  situations,  would  be  to  calculate  individual  stress  components 
explicitly  according  to  the  turbulent  Navier  Stokes  equation  first  established 
by  Osborne  Reynolds  in  1894.  The  solution  of  these  types  of  equations 
containing  turbulent  stress  terms  constitutes  the  multi-equation, 
stress-component  models.  This  type  of  turbulence  model  will  be  more 
universal  than  the  eddy-viscosity  type  of  models  but  requires  extensive 
computation.  This  will  be  the  model  of  the  next  decade. 

Formulated  using  mean-flow  velocity  components  u,  v,  w,  and  the  mean 
pressure  p,  the  Navier  Stokes  equation  modified  by  Reynolds  for  turbulent 
flow  resembles  the  original  equation.  The  right  hand  side  of  the  equation 
(when  n,  p  are  constants)  becomes: 
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In  the  formula,  the  instantaneous  velocity  component  (e.g.,  u)  is  the  sum 
of  the  mean  velocity  and  the  turbulent  velocity,  namely,  u  =  u  +  u’  in  the 
x-<iirection.  These  turbulent  velocity  components  (u*,  v’,  w’)  have  the 
property  that  their  long-term  (ensemble)  mean  value  approaches  zero.  But 
the  correlation  terms  such  as  u’v’,  etc.,  do  not  disappear  even  when  time 
averages  are  taken.  Reynolds  defined  these  terms  as  the  turbulent  stress 
terms.  They  are  as  follows; 


^xx  ^  ^  ’  ^yy  ^  ^  ’  ’'zz  ^  ^ 


The  Navier  Stokes  equation  modified  for  turbulent  flow  containing  the 
"Reynold  stress"  terms  looks  quite  similar  to  its  original  form: 
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(5.11) 


The  theory  of  turbulence,  as  well  as  the  turbulent  Navier  Stokes  equation 
in  its  original  form  proposed  by  Reynolds  in  1894,  has  several  difficulties  in 
application  because  of  the  following  (Hinze,  1959): 


•  When  the  partial  differential  equation  for  the  velocity  corrections  of  a 

given  order  are  derived  from  the  equations  of  turbulent  fluctuation, 
the  presence  of  the  inertia  terns  f-auses  the  appearance  of  the 
velocity  correlations  of  the  n.'..f,  .  order,  which  are  also 

unknown 

•  Equations  of  correlation  on  the  second  and  higher  orders  constructed 

out  of  the  equations  of  turbulent  fluctuation  contain  the  unknown 
terms  of  correlation  between  the  pressure  and  velocity  fluctuation 

•  The  value  of  the  decay  term  in  these  equations  has  to  be  determined. 


ALGEBRAIC  (ZERO-EQUATION)  MODEL  FOR  THIN-LAYER 
APPROXIMATIONS 

In  this  type  of  turbulence  closure  scheme,  the  eddy-viscosity  coefficients 
are  computed  based  on  the  Prandtl  formulation  in  which  the  mixing  length  is 
specified  algebraically.  Since  these  types  of  models  use  only  the  partial 
differential  equations  for  the  mean  field  and  no  differential  equation  for  the 
turbulent  quantities,  it  is,  therefore,  also  called  the  "zero  equation  model." 
These  models  relate  the  turbulent  shear  stress  only  to  the  mean  flow 
conditions  at  each  point  through  an  algebraic  relationship.  The  scheme 
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proposed  by  Baldwin  and  Lomax  (1978)  belongs  to  this  type  and  has  received 
the  most  attention  in  the  hypersonic  aerodynamic  flow  computation.  CFD 
models  for  simulating  the  SCRAMJET  flow  field  are  based  nearly  entirely 
(twelve  out  of  sixteen)  on  the  Baldwin-Lomax  scheme  (see  Table  4).  This 
scheme  has  several  major  advantages  in  computing  hypersonic  air  flow.  We 
will  discuss  this  closure  scheme  in  some  detail. 

The  basic  approach  of  the  Baldwin-Lomax  scheme  is  patterned  after  the 
Cebed  method  (1974).  However,  the  Baldwin-Lomax  scheme  has  the 
advantage  of  avoiding  the  necessity  for  finding  the  edge  of  the  boundary  layer. 
The  two-layer  model  divides  the  eddy  viscosity  into  two  parts,  as  follows: 


■  (/i^)  inner 
(/i^)outer 


y  <  y 

^  -  •'crossover 

y  <  y 

•'crossover  •' 


(5.12) 


where  y  is  the  local  normal  distance  from  the  wall  and  y-„..-  __  is  the 

Cl  OdoOVCA 

smallest  value  of  y  at  which  values  from  the  inner  and  outer  formulas  are 
equal.  In  other  words,  the  crossover  point  in  y  is  the  point  at  which 

becomes  less  than  (Mj.)jjjjjgj.-  Unlike  the  Prandtl  formulation  in  which  the 

eddy  viscosity  is  a  function  of  the  velocity  gradient,  the  Prandtl-Van  Driest 
formula  is  used  for  defining  the  eddy  viscosity.  In  that  formulation,  the 
vorticity  of  the  local  flow  field  is  used  in  the  inner  region: 


(u).  =  pP\ 0/1 

''^t'mner  ^  '  ' 


(5.13) 


The  mixing  length 


I  =  ky 


1-exp  ( 


-y+/A+) 


(5.14) 


where  y  is  the  normal  distance  from  the  wall,  y"*”  is  the  law-of-the-wall 
distance  (i.e.,  Vp  r  y/u  ),  k  is  the  von  Karman  constant,  and  is  a 

closure  contant  (A'^=  26).  The  magnitude  of  vorticity  |  ui\  is  defined  as  (for 
three-dimensional  flow''- 
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Table  4 

Use  of  Turbulence  Models  in  Recent  Aerospace 
Plane- Related  Hypersonic  Simulations 

Model  Numerical  Method  Turbulence  Closure  References 

1)  Supersonic  flow  2-D  explicit  scheme  2- layer  algebraic  eddy  Berman  and 

in  SCRAMJET  with  (MacCormack,  1969)  viscosity  (Baldwin  and  Anderson 

step,  H^injection  Lomax,  1978)  (1983) 

2)  2-D  analysis  2-D  explicit  scheme  2- layer  algebraic  eddy  Kumar 

of  SCRAMJET  in-  (MacCormack,  1969)  viscosity  (Baldwin  and  (1982) 
let  flow  field  Lomax,  1978) 

3)  Staged  in-  Time- split  algorit.  2- layer  algebraic  eddy  Veidner  A 

jection  for  (MacCormack  and  viscosity  (Baldwin  and  Drummond 

SCRAMJET  Baldwin,  1975)  Lomax,  1978)  (1982) 

4)  Hypersonic  lamin.  Predictor- corrector  2- layer  algebraic  eddy  Lawrence 

flow,  15®corner  implicit  (MacCormack,  viscosity  (Baldwin  and  et  al. 

at  M^  =  14.1  1982)  Lomax,  1978)  (1987) 

5)  NASA  modular  Time  unsplit/split  Baldwin  and  Lomax  White 

SCRAMJET  com-  MacCormack  and  Baldwin  (1978)  et  al.  (1987) 

bustor  flow  (1969,  1976) 

2-D,  PNS 

6)  3-D  hypersonic  Beam  and  Warming  Baldwin  and  Lomax  Thomas 

equilibrium  flow  ADI  (1976)  (1978)  (1988) 

with  ablation 


-64- 


Table  4  (Cont.) 

Use  of  Turbulence  Models  in  Recent  Aerospace 
Plane- Related  Hypersonic  Simulations 


Model 

Numerical  Method 

Turbulence  Closure 

Reference 

7)  Low- density 
hypersonic 
real  gas 

MacCormack 

implicit 

(1982) 

Baldwin  and  Lomax 
(1978) 

Hoffman  et  al. 
(1988) 

8)  Separated  flow, 
backflow 

turbulence 

Implicit  upwind- 
biased,  TVD,  appr. 

factorization 

Modified  Baldwin 
and  Lomax  (1978) 

Goldberg  and 
Chakravarthy 
(1988) 

9)  Separated  flow, 
backflow 

turbulence 

Implicit  upwind- 
biased,  TVD,  appr. 
factorization 

hybrid  k-L  one- 
equation  model 

Goldberg  and 
Chakravarthy 
(1989) 

10)  SCRAMJET  com¬ 
bustor  A  nozzle 
real  gas 

Implicit  finite- 
volume,  time¬ 
marching 

Baldwin  and  Lomax 
(1978) 

Nelson 

et  al.  (1989) 

11)  SCRAMJET  flow 
PNS,  combustor/ 
nozzle 

Implicit  upwind 

low  Re  k-  c 

two- equation 

Dash  et 
al.  (1989) 

12)  NASA  Ames  all¬ 
body  hypersonic 

Implicit  upwind 

PNS  solver 

Baldwin  and 

Lomax  (1978) 

Lockman  et 
al.  (1989) 

aircraft  model 
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I  >1  [  ^  j. 

I“l  =  l‘5r“S5rJ  +  [ 


dv  dw  V  ,  f  dw  5u  1 

"Si - 3FJ  l“3=i - 


(5,16) 


For  two-dimensional  flow, 


\u}\  -  \  ^  ^  1 

'  '  "  — srj 


(5.16) 


Eddy  viscosity  of  the  outer  region  is  a  function  of  the  Clauser  constant 
(C^p),  the  Klebanoff  intermittency  factor,  and  other  closure  constants: 


(^t^outer  ^  ^cp  ^  ^wake  ^Kleb 


(5.17) 


The  closure  constants  Baldwin  and  Lomax  used  have  been  determined  by 

requiring  agreement  with  the  original  Cebeci  formulation  for  constant 

pressure  boundary  layers  at  transonic  speeds.  They  used  C  =1.6, 

cp 

^Kleb~^'^’  K=0.0168.  is  defined  by  additional  parameters 

involving  the  difference  between  maximum  and  minimum  total  velocity  in  the 
velocity  profile.  is  usually  taken  to  be  the  minimum  of  ^max 

^wakeVix"’dif/''max  "^max  “  ^  ’'max 

with  F^^  being  the  maximum  value  of  F(y)  in  a  given  transverse  profile. 

The  function  of  F(y)  is  given  by 


F(y)  =  y|u|(l-exp(-y'''/A''')) 


(5.18) 


^diff  determined  by  taking  the  difference  between  the  maximum  and 
minimum  velocities  in  a  given  velocity  profile  with  min(U^£^)  equal  to  zero 
everywhere  else.  Therefore, 
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1 

Udiil=l"^  +  ''^l^  (519) 

The  value  of  represents  the  Klebanoff  intermittency  factor,  given 


by 


^Kleb^y)  “ 


1  -1-  5.5  [ 

^Kleb  y 

Y  J 

max 

-1 


(5.20) 


The  constants  for  the  Baldwin  and  Lomax  model  are;  A'^=  26.0,  C  = 

cp 

^kleb^  ^wake^  0.0168.  Prandtl  (Pr)  and 

turbulent  Prandtl  number  (Prt)  used  in  the  original  Baldwin-Lomax  model 
(1978)  are  0.72  and  0.9,  respectively. 

Since  the  model  uses  the  normal  distance  from  the  wall  as  its  primary 
parameter  in  determining  the  turbulent  viscosity.  If  certain  discontinuity 
(e.g.,  a  step)  exists,  there  will  be  a  jump  in  the  value  of  viscosity.  This 
ambiguity  sometimes  requires  smoothing.  This  can  be  done  by  a  relaxation 
equation  (Waskiewicz  et  al.,  1980,  Berman  and  Anderson,  1983)  of  the  form 


Mt-M, 


^e-  ^te 


te 


X  >  X 


te 


(5.21) 


where  subscript  t  represents  turbulent,  te  represents  trailing  edge,  le  denotes 
leading  edge,  6  is  the  boundary^ayer  thickness,  and  Ax  is  grid  spacing  of  the 
model. 

One  major  advantage  of  the  zero-equation  algebraic  turbulence  model  is 
the  simplicity  of  its  application.  It  can  be  modified  for  different  flow 
situations.  An  example  is  Goldberg’s  treatment  for  the  separated  flow  region 
(Goldberg,  1986,  Goldberg  and  Chakravarthy,  1988).  The  algebraic  backflow 
model  has  been  incorporated  into  a  Reynolds-averaged  Navier  Stokes  solver 
that  uses  the  Baldwin-Lomax  turbulence  model  outside  of  separation  bubbles. 
The  scheme  has  been  applied  to  calculate  the  reattaching  flow  over  a 
backward-facing  bump.  Data  comparison  indicates  that  this  combination 
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outperforms  the  more  complicated  e-t  (or  some  times  called  k-t) 
two-equation  model,  which  is  discussed  below. 


ONE-EQUATION  TURBULENCE  CLOSURE  MODELS 

In  zero-equation  models,  the  closure  is  made  through  the  mean-velocity 
field.  The  specification  of  turbulence  field  by  means  of  an  algebraic 
relationship  implies  that  generation  and  dissipation  of  turbulent  energy  are  in 
balance  everywhere.  Therefore,  the  dynamic  process  of  convection  and 
diffusion  of  turbulent  energy  is  ignored.  These  are,  of  course,  important 
features  of  real  turbulent  flows.  To  include  these  processes,  the  transport  of 
turbulent  energy  has  to  be  computed  explicitly  via  another  set  of  partial 
differential  equations.  In  fact,  at  the  very  beginning,  the  founders  of  modern 
tubulence  theory  (namely,  Prandtl  and  Kolmogorov)  suggested  that  turbulent 
viscosity  be  determined  by  way  of  differential  rather  than  algebraic  equations. 

In  the  set  of  differential  equations  the  dependent  variable  often  being 
computed  for  the  closure  is  the  sub-grid-scale  (SGS)  turbulent  energy  density 
"e"  (or  sometimes  called  "k")  within  a  computational  cell.  In  most  models, 
the  following  transport  equation  for  e  is  used: 


de 

dx 

(1) 


1  TT  —  9  ^t  de 
'^i  dx..  ~  dx.-  CT  dx- 

1  ^  L  ®  * 

(2)  (3) 


5U: 


(4) 


(5) 


(5.22) 


where 

(1)  =  local  acceleration; 

(2)  =  convection  of  sub-grid-scale  turbulent  energy; 

(3)  =  diffusive  transport  of  SGS  energy; 

(4)  =  energy  production  by  shear  stress; 

(5)  =  viscous  dissipation; 

=  turbulent  viscosity,  which  is  assumed  to  be 

a  property  of  the  local  state  of  the  turbulence; 
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=  effective  Prandtl  number  for  the  diffusion  of 
e 

turbulent  energy. 

In  most  one-equation  models,  the  turbulent  viscosity  is  determined  as  a 
function  of  the  local  energy  intensity  and  the  local  mixing  length  (/)  from  the 
boundary  via  the  following  relationship 

Mt  =  C^/^/Re  (5.23) 

The  rate  of  viscous  dissipation  is  determined  by  the  local  energy  density 
according  to  the  following  relationship, 

3 

f  =  Re  (5.24) 

where  is  one  of  the  closure  constants  (approximately  0.9),  and  /  is  the  local 

mixing  length,  which  is  usually  proportional  to  the  distance  from  the  solid 
wall. 

The  sequence  of  closure  computation  also  involves  the  relationship 
between  pressure,  density,  and  compressibility,  which  is  approximated  by  the 
equation  of  state  of  air.  To  explain  the  cyclic  process  of  numerical  integration 
in  which  "one"  additional  partial  differential  equation  is  used  for  SGS  energy, 
the  sequence  starting  from  an  arbitrary  point  in  time  during  the  process  goes 
as  follows. 

The  velocity  fields  are  computed  in  the  entire  domain  by  the  balance  of 
pressure  gradient,  shear  stress,  and  other  forces  including  the  specified 
far-field  boundary  conditions.  From  these  velocity  fields  just  determined, 
heat,  SGS  energy,  and  other  essential  constituents  are  transported  and 
diffused.  From  the  transport  of  SGS  energy,  the  new  turbulent  viscosity 
coefficients  throughout  the  entire  flow  field  are  determined.  From  these  stress 
coefficients  and  other  diffusion  coefficients  computed  in  a  similar  manner, 
together  with  the  updated  pressure  field,  a  new  velocity  field  is  computed. 
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From  the  computed  velocity  and  temperature  field,  the  drag,  lift,  detailed 
skin  friction,  and  heat  coefficients  about  the  entire  vehicle  can  be  computed. 


TWO-EQUATION  TURBULENCE  CLOSURE  MODELS 


In  essence,  the  one-equation  turbulence  model  is  based  on  the  proposal 
made  earlier  that  the  sub-grid-scale,  time-averaged,  turbulent  kinetic  energy 
e  (or  sometimes  called  k) 


+  V 


,2 


+  w 


,2 


(5.25) 


is  determined  from  the  solution  of  a  convective  transport  partial  differential 
equation  similar  to  the  transport  of  heat.  In  the  above  equation,  u’,  v’  and  w’ 
are  turbulent  "random  fluctuations"  around  the  time-averaged  mean 
velocities.  The  major  shortcomings  of  the  one-equation  model  are:  1)  the 
influence  of  convection  and  diffusion  on  the  small-scale  turbulent  random 
velocities  is  not  accounted  for;  and  2)  the  viscosity  vanishes  whenever 

du/dy  is  zero.  To  avoid  these  problems  we  need  one  more  equation  to 
describe  the  transport  of  the  length  scale  1.  The  turbulent  viscosity  is 

determined  in  the  same  way  as  in  the  one-equation  model,  namely, 

1 

Mt  =  ye^l  (5.26) 


In  other  words,  the  length  scale  /  is  calculated  by  a  differential  equation 
rather  than  prescribing  it  algebraically.  In  most  two-equation  models,  the 
length  scale  is  determined  indirectly  via  a  variable  z 


m  fl. 
z  =  e  r 


(5.27) 
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where  m  and  n  are  constants  that  vary  according  to  the  model  proposed  by  a 
different  modeler,  as  shown  in  Table  5. 


Table  5 

Closure  Constants  Proposed  by  Different  Modelers 
in  Two-Equation  Turbulence  Modeling 


Proposer  z 

m,  n 

Kolmogorov  (1942)  z  =  e  1 

m  =  — 2“  ,  n  =  — 1 

3 

Chou  (1945)  z  =  e^  1 

m  —  — 2“  ,  n  —  —1 

Rotta(1951)  z  =  l 

Spalding  (1967a) 

3 

II 

o 

B 

II 

Rodi  and  Spalding  (1970)  z  =  e  1 

m  =  1,  n  =  1 

Spalding  (1969)  z  =  e 

3 

II 
»— • 

B 

II 

The  transport  equation  of  e  and  z  can  be  rearranged  to  have  the  similar 
form 


(5.28) 


f  ''t  “  1  j. . 

>t 

r  ^n2 

e 
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Dz 

d 

["t 

dz] 

+  z 

[c 

1 

T>t““ 

ik. 

ik.\ 

*^£1  e 

. 

(5.29) 


In  these  equations,  c^,  cr^,  c^j,  c^2>  6*^.,  are  turbulent  closure 

constants.  They  have  to  be  determined  by  fitting  them  against  experimental 
data.  A  set  of  optimally  fitted  turbulent  closure  constants  are  listed  in  Table 
6. 


Conventional  two-equation  modds  formulated  using  e-z  (or  or  sometimes 
called  k-€  )  generally  are  inaccurate  for  boundary  layers  in  an  adverse 
pressure  gradient.  Wilcox  has  recently  (1988)  proved  that  the  use  of  "wall 
functions"  (Rodi,  1981)  tends  to  mask  the  shortcomings  of  such  models  and  is 
inadequate  for  flows  with  mass  injection.  Under  the  latter  case,  skin  friction 
predicted  by  the  k-€  modd  is  as  much  as  50  percent  higher  than  measured 
(Wilcox,  1988).  Using  a  singular  perturbation  method,  Wilcox  proposed  a 
multiscale  model  with  no  viscous  damping  of  the  model’s  closure  coeffidents 
which  gives  improved  accuracy  for  flow  with  surface  mass  addition. 

Strictly  speaking,  the  one-equation  modeling  approach  is  not  as  general  as 
the  two-equation  modd.  In  practice,  however,  this  method  allows  a  modder 
to  select  an  algebraic  function  for  the  length  scale  according  to  the  specific 
physical  process  guided  by  experimental  data  and  field  observations.  The 
one-equation  approach  thus  avoids  the  problematic  z-equation  (see  Goldberg 
and  Chakravarthy,  1989),  which  is  needed  to  prescribe  the  length  scale  in  the 
two-equation  model.  This  advantage  was  reported  by  Liu  and  Leendertse 
(1978,  1987,  1990)  in  three-dimensional  CFD  Navier  Stokes  models  of 
geophysical  fluid  dynamic  systems. 
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Table  6 

Optimally  Fitted  Turbulence  Closure 
Constants  in  CFD  Modeling 


c 

fi 

^el 

^e2 

Reference 

0.09 

1.00 

1.30 

1.44 

1.92 

Launder  R  Spalding  (1974) 

0.09 

1.00 

1.30 

1.44 

1.92 

Launder  k  Sharma  (1974) 

0.09 

1.00 

1.30 

1.45 

2.00 

Hassid  and  Poreh  (1978) 

0.09 

2.00 

3.00 

1.81 

2.00 

Hoffman  (1975)^ 

0.09 

0.90 

0.95 

1.35 

2.00 

Dutoya  t  Michard  (1981) 

0,09 

1.00 

1.30 

1.35 

1.80 

Chien  (1982) 

0.084 

1.69 

1.30 

1.00 

1.83 

Reynolds  (1976)^ 

0.09 

1.00 

1.30 

1.43 

1.94 

Launder  et  al.  (1973) 

0.09 

2.00 

(Not  directly  compatible) 

Vilcox  and  Rubesin  (1980) 

0.09 

1.00 

1.30 

1.44 

1.92 

Warfield  and  Lakshminarayana 

(1987) 

a 


Model  constants  have  to  satisfy  the  relationship 
2  1/2 

c^j  =  c^2  ■  ^  ^  equation  reduces  to  zero. 
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STRESS-COMPONENT  CLOSURE  MODELS 


At  the  present  time,  the  most  complete  set  of  turbulent  closure  equations 
is  the  "mean  Reynolds  stress"  approach,  which  was  proposed  by  Chou  (1945) 

and  Davidov  (1959,  1961).  The  formulation  includes  the  transport  and 

_ _ 2 

diffusion  of  length  scale  1,  u’v’,  u’^,  u’v’  terms,  and  results  in  3  equations  for 
mean  flow  components,  1  equation  for  continuity,  3  equations  for  shear 
stresses,  3  equations  for  normal  stresses,  10  equations  for  triple  correlations, 
for  a  total  of  20  equations  for  a  homogeneous  aerodynamic  systems.  The 
turbulent  stress  component  formulated  originally  by  Chou,  in  compact  tensor 
notation,  for  double  velocity  correlation,  takes  the  following  form: 


p  5t  p 


P 


P 


Pi’^k 


+  P 


,k''i 


V 

P 


V^T- 


ik 


U: 


i,m'^k,n. 


(5.30) 


Ten  equations  for  triple  correlations  are: 


^iVl 


k - +  Uj  j  uJu-Uj  -H  Uj^  j  uJujUj  -h  Ui  j  u\uj^  -h  Uj  uJujU^  -h 


U^h\“lJ,j+[^SVlj,j  =  -  - f  [p,iVl  +  P,kVi  +  P,lVk 


I  U  Ui  U,  I 

[  1  k  IJ  ,mn 


u.  u,  Ui -f  Ui  Ui  u- -I- u,  u-  Ui 
°  (  i,m  k,n  1  k,m  l,n  i  l,m  i,n  li 


(5.31) 


The  continuity  equation  is: 
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u^.  =  0 
)J 


(5.32) 


Since,  in  tensor  notation,  a  subscript  preceded  by  a  comma  denotes  the 
covariant  derivative,  the  continuity  equation  in  scaler  notation  is  then. 


^1  I  ^2  , 

£hci  dni 


(5.33) 


defining  q^  =  UjU^ 


=  the  variance  or  the  RMS  of  the  velocity  fluctuation, 


7J(x‘xV),  (j=l,2,3). 

The  triple  and  quadruple  correlation  technique  is  essentially  a  method  of 
successive  approximation  to  the  solution  of  the  turbulence  problem  by  solving 
(PDEs)  for  each  component  of  Reynolds  stress.  The  approximation  takes  the 
following  order; 


Initial  approximation;  Reynolds  equations  of  mean  motion,  which  contain 
the  unknown  apparent  stresses 

Second  approximation;  Solve  the  equations  of  mean  motion  and  of  the 
double  correlation  by  making  certain  approxima¬ 
tions  to  the  triple-velocity  correlations  in  the 
equations 

Third  approximation;  Solve  the  equations  of  mean  motion  and  both  the 
double  and  triple  correlations  simultaneously  by 
assuming  approximations  for  the  quadruple  correla¬ 
tions. 


The  above  implies  that  there  will  always  be  more  unknowns  than  the 
available  equations.  Or  one  can  obtain  an  approximation  at  a  certain  level 
with  some  remaining  lower-order  unknown  terms,  thus  the  term  "turbulence 
closure." 


For  three-dimensional  turbulence  closure  computations,  in  addition  to  the 
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dynamic  and  chemical  balance  equations,  there  are  six  components  of  the 
Reynolds  stresses  u^Uj  and  three  flux  (scalar)  components  plus  one  transport 

equation  for  the  scalar  concentration  fluctuations,  for  a  total  of  10  additional 
partial  differential  equations  to  solve.  The  application  of  the 
stress-component  turbulence  modeling  is,  therefore,  still  very  limited  due  to 
its  computational  demand.  Some  simplification  can  be  made  to  the  PDE  such 
that  they  reduce  to  algebraic  expressions  but  still  retain  the  fundamental 
characteristics  of  the  basic  approach.  In  the  transport  equations,  if  the 
gradients  of  the  dependent  variables  (e.g.,  rate  of  change,  convection,  and 
diffusion)  are  eliminated  by  model  approximations,  the  differential  equations 
can  then  be  converted  into  algebraic  expressions  (Launder  et  al.,  1975; 

Launder  and  Spaulding,  1974).  Warfield  and  Lakshminarayana  (1987)  made 
some  tests  with  this  approach. 

Difficulties  associated  with  the  length-scale  equation  often  cause  the  poor 
performance  of  the  two-oquation  turbulence  models.  Since  the  stress 
component  model  uses  a  similar  length-scale  equation,  this  type  of  model  may 
not  necessarily  give  better  predictions  over  the  one-  or  two-equation  model. 
However,  this  remains  to  be  seen. 


TWO-FLUID  MODEL  OF  TURBULENCE  AND  COMBUSTION 

In  many  turbulent  combustion  processes  (e.g.,  within  a  SCRAM  JET 
engine),  two  interacting  fluids  of  different  states  (e.g.,  velocities)  are  assumed 
to  exist  within  a  finite-control  volume.  The  dynamic  and  chemical  processes 
of  the  interspersed  gas  fragment  within  the  flame  can  sometime  be 
represented  more  faithfully  by  using  a  two-fluid  model  of  the  flame  than  a 
single  fluid  formulation. 

Introduced  originally  by  Shchelkin  (1943)  and  Wohlenberg  (1953),  the 
theory  of  two-fluid  model  of  turbulence  has  recently  been  advanced  to  use  the 
analytical  tools  for  two-phase  flows  (Spalding,  1986;  Fan,  1988).  If  R 
represents  the  volume  fraction,  the  mean  density  of  the  fluid  is 
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P  ~  ^2^2  (5.34) 

When  analyzing  the  two-fluid  flow,  the  averaged  fluid  property  ^  can  be 
defined  in  two  ways: 

1.  the  time  (or  volume)  average:  <})  =  Rj(|)j  -I-  R2<j)2  (5.35) 

2.  the  mass-weighted  average:  (j)  =  (RjPj(j>j  -h  R^p2^2)  (5.36) 

J 

The  momentum  transfer  or  interfluid  friction  in  the  turbulent  closure  can 
be  expressed  as  fj2  : 

=  C(  a-'/  I  Vj-Vj  I  (6.37) 

where  f ^2  represents  the  interfadal  friction  between  the  two  fluids,  c^  is  a 

dimensionless  proportionality  constant,  a~^  is  the  amount  of  fluid/fragment 
interface  area  per  unit  volume  of  space,  p  is  the  density  of  the  lighter  of  the 
two  fluids,  and  I VJ-V2I  is  the  local  time-averaged  relative  speed  of  the  two 

fluids.  The  area/volume  quantity  (i.e.,  a~^)  is  taken  as  being  proportional  to 
the  volume  fraction  product  RjR2>  so  that  it  vanishes  when  either  fluid 

disappears.  During  the  dissipation  of  turbulence,  the  turbulent  fluid  can 
transfer  to  the  nonturbulent  fluid.  On  the  other  hand,  the  nonturbulent  fluid 
can  become  turbulent  as  a  result  of  the  volumetric  entrainment  process. 

The  basic  governing  equation  for  the  mean  flow  contains  two  fluids 
co-existing  in  one  space; 

VA) .  «  f  „  .  «fk  1 

ac.  "'iT 

J  J  J  ■* 
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dx.- 
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A  D.,  - - 
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+  V  +  ’i).k 


(5.38) 


where 

k  =  subscript  denoting  fluid  1  or  fluid  2 
j  =  subscript  representing  coordinate 

r  =  turbulent  or  laminar  diffusion  coefficient  within  one  fluid 
D  =  coefficient  of  diffusion  transport  of  dependent  variable 
due  to  the  fluctuations  of  velocity  and  volume  fraction 
S  =  within-fluid  source  term 
I  =  inter-fluid  source  term 

At  present,  the  fluid  model  has  been  formulated  for  analyzing  a  tubulent 
premixed  combustion  process.  Documented  computer  programs  for  one-  and 
two-phase  flows  such  as  the  PHOENICS  code  are  available  from  the  London 
Imperial  College. 


MODELING  COMPRESSIBLE  TURBULENCE 

In  modeling  turbulence,  the  compressibility  effect  can  generally  be 
neglected  if  the  ratio  between  time-averaged  air  density  and  the  variation 
of  density  associated  with  the  turbulent  fluctuation  (p)  is  small,  i.e., 

p  /  p  <<  1  (5.39) 

This  is  approximately  the  same  order  of  u’7  c*  (Hinze,  1959),  which  is 
the  square  of  the  Mach  number  of  the  turbulence.  In  compressible  turbulent 
flow,  the  mean  motion  is  affected  not  only  by  turbulence  shear  stresses  but 
also  by  stresses  generated  from  double  and  triple  correlations  involving  the 
density  fluctuation  p.  In  hypersonic  compressible-flow  jets,  variation  in 
density  occurs  from  pressure  differences.  Only  limited  information  is 
available  on  the  effect  of  compressibility  of  air  on  the  mechanism  of  turbulent 
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air  flow  in  a  boundary  layer.  Near  boundaries,  the  reduction  in  high  velocity 
causes  a  conversion  into  heat  through  compression.  Heat  generated  by  the 
compressional  process  does  not  distribute  uniformly  within  the  boundary 
layer.  The  temperature  variation  due  to  this  process  changes  fluid’s 
properties.  A  transport  of  heat  occurs  by  the  molecular  and  turbulent 
diffusion  processes.  The  turbulence  eddy  entrainment  rate  of  the  compressible 
turbulent  flow  is  less  than  that  of  the  corresponding  incompressible  flow.  The 
effect  of  the  shock  wave  on  the  rate  of  energy  transport  into  or  out  of 
turbulence  eddies  is  not  well  understood  at  this  time.  For  compressible 
turbulence,  the  fundamental  spectral  laws  that  govern  the  partition  and  the 
cascade  of  turbulent  energy  as  proposed  by  Kolmogorov  may  need  revision 
(Dimotakis,  1989).  Recent  research  on  compressible  turbulence  includes 
mainly  theoretical  studies  such  as  the  three-dimensional  structures  of 
compressible  shear  layer  (Papamoschou,  1989),  and  the  theory  of  sonic  eddy 
(Breidenthal,  1990).  In  the  latter  theory,  the  effect  of  Mach  number  on 
turbulent  shear  flow  is  expressed  via  the  concept  of  sonic  eddy.  Thus, 
obviously,  it  may  be  some  time  before  our  understanding  of  compressible 
turbulence  models  can  be  refined  and  verified,  and  proper  models  used  to 
calculate  hypersonic  flows.  This  will  require  a  combination  of  measurements, 
analyses,  and  theory. 


AREAS  OF  LOW  PREDICTIVE  RELIABILITY  IN  MODELING 
TURBULENCE 

In  this  section,  we  have  reviewed  the  necessity  for  modeling  turbulence  in 
the  CFD  simulation  of  hypersonic  aerodynamic  flows.  In  some  cases,  the 
basic  theory  of  turbulent  flow  per  se  caimot  truthfully  describe  the 
complicated  physical  phenomena  of  fluid  turbulence.  In  some  cases,  we  are 
limited  by  the  capacity  of  the  present  day  computer,  so  that  CFD  models  do 
not  have  high  enough  resolution  to  simulate  turbulent  flow  at  high  Reynolds 
number.  Other  difficulties  involve  the  lack  of  sufficient  data  to  verify  many 
turbulent  models  such  as  those  proposed  just  for  the  simpler  two-equation 
model  (as  tabulated  in  Table  5).  As  we  have  quoted  at  the  beginning  of  this 
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section,  the  NRC  conclusions  on  the  present  status  of  turbulence  modeling 
suggest  that  extensive  research  is  still  needed  before  we  have  sufficient 
confidence  in  CFD  predictions  particularly  in  areas  without  verification  data. 

From  our  study,  we  have  summarized  a  list  of  difficult  subjects  and  areas 
with  low  predictive  reliability  in  turbulence  flow.  Table  7  lists  these.  Several 
of  these  areas  are  associated  with  anisotropic,  nonhomogeneous  turbulence, 
which  induces  hydrodynamic  instability.  The  combustion  process  in  a 
SCRAMJET  engine  may  involve  nonhomogeneous  tubulence  if  steps  or  blocks 
are  used  for  inducing  more  complete  mixing  at  such  high  speed.  The  problem 
of  hydrodynamic  stability  has  been  studied  for  decades  by  many  well-known 
scientists  in  the  fluid  dynamic  field.  We  do  not  yet  have  a  final  conclusion  on 
the  basic  stability  criterion  for  the  nonhomogeneous  turbulence  problem  (see 
Table  8).  It  has  long  been  recognized  that  energy  dissipation  in  turbulence 
must  be  intermittent  in  space.  Recently,  the  spectral  method  has  been  used 
in  conjunction  with  supercomputers  to  study  the  intermittency  in  turbulence 
with  2  million  (128^)  grid  resolution  (Hosokawa  and  Yamamoto,  1990). 

Direct  simulation  methods  have  also  been  used  to  study  the  effect  of  Mach 
number  on  the  stability  of  plane  supersonic  wake  (Chen  et  al.,  1990),  so  that 
the  physics  of  linear,  nonlinear,  and  three-dimensional  stages  of 
laminar-turbulent  free-shear  flow  transition  is  better  understood.  In  the  area 
of  mixing  and  combustion,  in  incompressible  flows,  xhe  mixing  layer  is 
generally  convectively  imstable.  Similar  information  is  not  available  for  the 
compressible  mixing  layer  for  which  temperature  effects  are  important. 

Recent  experimental  results  (Jackson  and  Grosch,  1990)  suggest  a  way  to 
derive  from  linear  stability  theory  a  "convective  Mach  number"  (proposed 
originaly  by  Bogdanoff,  1983)  for  a  compressible  mixing  layer  for  multispecies 
gas.  This  provides  a  convenient  way  for  studying  the  stability  of  compressible 
free  shear  layer.  In  supersonic  combustion  jet  engine  design,  it  is  important 
to  understand  the  behavior  of  compressible  turbulent  shear  layers.  With 
progress  in  numerical  modeling,  computing  capability,  and  more  experimental 
data,  our  understanding  of  fluid  turbulence  will  certainly  improve,  which 
would  also  lead  to  higher  predictive  reliability  for  supporting  aircraft  design. 

A  key  unresolved  issue  is  the  analysis  of  the  laminar-turbulent  transition 
process,  first  for  the  low-speed  case  and  ultimately  for  the  high-speed  case  of 
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Table  7 

Modeling  Turbulence: 

Difficult  Subjects  and  Areas  Vith  Low  Predictive  Reliability 


Strong  aerodynamic  curvature 

Intermittency  and  large-  scale  flow  structure 

Rapid  compression- expansion  (k- £  model  needs  modifications) 

Kinematically  influenced  chemical  reaction 

low  Reynolds  number  effects 

Strong  swirl 

Aerodynamic  turbulence  is  strongly  influenced  by  body  force  acting 
in  a  preferred  direction  (k- c  model  no  longer  valid) 
Uncertainties  in  setting  the  boundary  condition 
Large- density  fluctuations  ->  high  Mach  number 
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Table  8 


Modeling  Nonhomogeneous  Turbulence: 
Stability  Criterion  for  Stratified  Shear  Flow 


»cr 

Proposed  by 

Year 

Conditions  and  Assumptions 

1.0 

force) 

Richardson 

1920 

(Roy.  Soc.  20,354) 

Boussinesq  approximation  (deviation 
of  p  is  neglected  except  in  buoyancy 

Eddy  diffusion  for  heat=momentum 

1.0 

Taylor 

1931  (1915) 

(Sci.  Paper  2,240) 

Same 

2.0 

Prandtl 

1930 

(Springer) 

Taylor  in  1931  suggested  multiplying 
by  0.5  to  obtain  1.0 

1/4 

Miles 

1961 

(JFM  10,496) 

Dynamically  attainable 

motions  using  energy  consideration 

1/4 

Chandrasekhar 

1961 

(Clarendon,  Oxford) 

Same 

1/4 

Howard 

1961 

(JFM  10,  509) 

Same 

1.0 

Abarbanel,  Holm 
and  Ratiu 

,  1986 

(Roy.  Soc.  318,349) 

Proved  R  >  1  is  sufficient 
for  Liapunov  stability  of 

3-D  stratified  shear  flow 

Boussinesq  approximation 

1.0 

Miles 

1986 

(Phy.  Fid.  29,3470). 

Amended  Chandrosekhar’s  flaw 
in  energy  derivations  of  1961 

* 

Eichardson  number 


^  Critical 
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interest  to  NASP  designers.  At  present,  there  is  little  analysis  or  data  to 

guide  o.  prediction  of  transition,  except  for  the  use  of  linear  stability  and  the 
N 

classic  e  method  that  has  been  useful  in  lower-speed  flows.  Unfortunately 
the  e^  method  gives  little  useful  detail  beyond  the  location  of  transition,  and 
presumably  only  for  flows  that  have  disturbance  similar  to  those  that  were  in 
the  original  database  for  determining  "N".  Many  analysts  have  suggested 
that  uncertainties  in  the  location  of  transition  on  the  forebody  of  a  NASP 
vehicle  could  entail  uncertainties  of  as  much  as  200  percent  in  vehicle  weight. 
What  is  also  suspected  is  that  the  spatial  and  temporal  characteristics  of  the 
transition  process,  in  terms  of  intermittency,  unsteadiness,  and 
three-dimensionality,  and  even  the  presence  of  a  large  unsteady  flow 
structure,  can  markedly  affect  inlet  behavior  and  the  ability  of  the 
forebody-inlet  combination  to  perform  properly.  Tr?'  -ition  is  still  an 
experimental  art  form  and  there  are  limited  da*  ;mde  analysts.  It  may  be 
years  before  we  obtain  the  needed  statistical  data  to  properly  determine  the 
characteristics  of  hypersonic  flow  fields.  Other  issues  involve  the  possibility 
of  relaminization  in  high-speed  flow  and  the  behavior  of  the  hypersonic 
turbulent  shear  layers. 

Density  fluctuation  should  play  an  important  role  in  very-high-speed 
turbulent  boundary  layers,  but  there  are  almost  no  useful  data  to  guide  the 
formulation  of  CFD  models  reflecting  such  phenomena.  In  addition,  little  is 
known  about  the  statistics  of  turbulent  fluctuations  and  the  occurrence  of 
large-scale  fluid  structures  in  hypersonic  flow.  Such  large  eddy  structures 
exist  in  other  types  of  turbulent  boundary  flow.  Even  more  troublesome  than 
our  present  lack  of  understanding  of  the  hypersonic  turbulent  boundary  layer 
is  the  lack  of  data  and  understanding  concerning  the  hypersonic  turbulent 
mixing  process  that  occurs  in  the  engine.  Current  data  are  limited  but 
suggest  that  the  efficiency  of  mixing  of  two  streams  decreases  rapidly  as  the 
Mach  number  of  the  high-speed  stream  increases,  that  large-scale  fluctuations 
are  likely  to  occur,  and  that  compressible  shear  layers  may  be  far  more  stable 
than  previously  thought.  Thus,  shear  layers  may,  if  left  to  their  devices,  be 
more  laminar  than  is  desirable.  This  suggests  that  means  may  be  needed  to 
augment  compressible  mixing,  perhaps  using  small  eddy-generating  devices. 
These  devices  might  produce  small  shock  waves  (shocklets)  that  could  assist 
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eddy  formation  and  promote  mixing.  The  resulting  temporal  fluctuations  in 
the  engine  flow  path  could  influence  inlet  performance  by  feeding  disturbances 
upstream  through  the  thick  boundary  layers,  leading  to  inlet  buzz  or  unstart. 

As  far  as  inlet  performance  is  concerned,  there  is  also  little  evidence  that 
inlet  dynamics,  which  can  be  affected  by  the  forebody  flow  upstream  and  the 
combustor  flow  downstream,  can  yet  be  modeled  by  CFD  in  a  realistic  way. 
Furthermore,  the  need  to  operate  SCRAMJET  at  full-length  scale  to  allow 
for  complete  mixing,  and  the  effects  of  combustion  instabilities  and  other 
fluctuations  on  inlet  unstart,  suggest  that  without  a  full-scale  test  that 
includes  combustion,  we  will  not  have  confidence  in  our  ability  to  diagnose 
the  cause  of  inlet  malfunction.  All  of  these  performance-related  questions 
need  to  be  addressed  by  CFD  modeling  in  the  near  future. 
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6.  SUPERCOMPUTERS,  PARALLEL  PROCESSING,  AND  VECTOR 

PROGRAMMING 


INTRODUCTION 

Most  NASP-rdated  CFD  simiilations  use  supercomputers.  This  section 
attempts  to  evaluate  the  following  aspects  of  NASP-related  subjects  in  super 
computing: 

(1)  NASP-related  CFD  modeling  in  the  past  two  decades. 

(2)  Present  hardware/software  capability  and  costs. 

(3)  Technical  limitations  in  hardware,  software,  and  the  future  trend  of 
hardware/software  development  available  for  NASP  research/design 
support. 

Looking  back  in  the  history  of  computers,  numerical  weather  prediction 
and  nuclear  reactor  simulation  are  two  of  the  driving  forces  for  the 
development  of,  and  also  major  users  of,  the  early  computers.  Both  involve 
numerical  solutions  of  the  Navier  Stokes  equations.  At  each  stage  of  the 
development,  these  simulations  were  performed  on  the  largest  computers 
available  at  the  time.  However,  most  of  those  large  computers  were  "scalar" 
machines  designed  for  general  purposes,  particularly  for  business  data 
processing.  One  of  the  earliest  group  of  numerical  simulations  involving  the 
solution  of  the  "parabolized"  Navier  Stokes  equations  (PNS)  for  hypersonic 
aerodynamic  flows  were  conducted  at  RAND  using  general  purpose  IBM 
business  computers  (Cheng  et  al.,  1970).  Using  scalar  machines  for  solutions 
of  vector  quantities,  the  integration  process  is  carried  out  sequentially,  which 
is  not  as  efficient  as  if  it  were  conducted  in  parallel. 

Special  machines  were  designed  based  on  this  parallel  principle. 
ELLIAC-4,  which  consists  of  64  parallel  processors,  is  one  of  the  earliest 
examples.  But  these  special-purpose  machines  were  expensive  and 
time-consuming  to  build  since  they  were  not  mass-produced. 
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SUPERCOMPUTERS  AND  PARALLEL  PROCESSING 

Since  the  mid-TOs,  manufacturers  such  as  Control  Data  and  Texas 
Instruments  began  to  build  machines  for  scientific  computations.  Control 
Data  started  unth  the  STAR  series,  which  was  the  predecessor  of  its  Cyber 
series  vector  machines.  Since  the  introduction  of  CRAY-1  in  1976,  the  total 
number  of  "supercomputers"  has  grown  to  a  total  of  180  machine  from 
several  computer  makers.  Most  hydrodynamic  codes  run  on  these  machines, 
including  NASP-related  hypersonic  simulations. 

Hypersonic  simulations  run  faster  on  vector  machines.  Technically,  the 
major  reason  is  that  in  traditional  sequential  machines,  the  system  contains 
a  memory,  an  instruction  processor,  an  arithmetic  processors  and  an 
input/output  system.  In  the  memory,  data  and  instructions  occupy  unique 
addresses.  Operations  instructions  have  to  carry  addressing  information  to 
access  the  required  operands  in  the  memory.  The  final  computational  speed 
is  determined  by  the  slowest  component  that  controls  the  machine’s  cycle 
time.  Vector  processors,  on  the  other  hand,  operate  by  a  pipeline  procedure 
that  performs  many  arithmetic  operations  concurrently  with  a  single  vector 
instruction,  thus  achieving  much  higher  system  throughput. 

Data-swapping  between  CPU  and  the  external  memory  units  is  also  a 
routine  operation  in  CFD  modeling.  During  numerical  integration  of  the 
Navier  Stokes  equations,  the  portion  not  being  actively  handled  is  usually 
stored  temporarily  in  the  external  memory  unit.  Traditionally,  external 
magnetic  disk  units  are  used.  Some  supercomputers,  such  as  the  X-MP, 
offer  solid-state  storage  devices  as  an  alternative,  which  are  about  30  times 
faster  than  the  magnetic  units.  External  memory  units  are  usually 
configured  in  blocks  of  8-mega  bits.  The  efficiency  and  overall  throughtput 
of  a  CFD  simulation  depends  on  a  combination  of  hardware  and  the 
simulation  program  in  which  the  block  size  of  the  data  transfer  usually  also 
matters. 

In  February  1988,  Cray  introduced  the  model  Y-MP/832,  which  offers 
two  to  three  times  the  performance  of  X-MP.  The  Y-MP  has  8  CPUs  and  a 
central  memory  capacity  of  32  million  words.  With  eight  CPUs,  the 
Y-MP/832  has  twice  the  number  of  processors  and  memory  as  the  X-MP. 
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Each  CPU  operates  on  a  6-nansecond  clock  cycle.  At  clock  speeds  of  six 
billionths  of  a  second,  the  machine  can  thus  handle  32  million  64-bit  words 
in  central  memory.  Abour  ten  Y-MP  were  delivered  in  1989.  The  other 
supercomputer,  Control  Data’s  CYBER-205,  has  a  cycle  time  of  20  ns  which 
is  about  half  the  speed  of  Cray’s  X-MP. 

A  new  machine  (i.e.,  ETA-10)  introduced  in  Control  Data  is  comparable 
to  the  Cray’s  Y-MP.  The  ETA-10  was  introduced  in  June  1987.  The 
system  is  a  highly  vectoring,  parallel-processor,  with  virtual  memory.  A 
major  difference  between  ETA-10  and  other  supercomputers  is  the 
instruction  set  of  its  central  processing  unit.  Unlike  most  supercomputers 
which  increased  speed  by  reducing  the  instruction  set  (RISC),  the  ETA  10 
can  be  qualified  as  a  VCISC-a  very  complex  instruction  set  computer. 
Designed  mainly  for  large-scale  computational  fluid  dynamic  and  structure 
applications,  consequently,  the  use  of  vectors  is  inherent  both  in  the 
hardware  and  the  instruction  set.  Each  CPU  contains  in  addition  to  the 
scalar  processor  a  dual-pipe  line  vector  processor,  both  of  which  operate 
comoletely  in  parallel.  Each  CPU  is  a  44-layer  printed  single  circuit  board. 
The  ETA  10  system  can  contain  one  to  eight  processors  each  with  a  32 
megabytes  of  local  memory.  The  performance  benchmark  for  the 
single-processor  model-P  (air-cooled,  24-ns  clock)  is  equivalent  to  44  times 
the  speed  of  the  popular  Micro  VAX  II.  The  faster  model  ETA  10-G  is 
liquid  nitrogen  cooled  with  a  7-nanosecond  clock  and  has  a  rated  speed  of  94 
MFLOPS  when  equipped  with  a  single  CPU.  The  delivery  date  for  ETA 
10-G  was  in  December  1988.  The  price  basic  configuration  is  $1  M  for  the 
24-ns  ETAIO-P  and  $13.5  M  for  the  7-ns  ETAIO-G.  Because  of  financial 
difficulties,  the  production  of  ETA-10  machines  stopped  in  the  spring  of 
1990. 

The  next  model  being  developed  at  Cray  Research  is  said  to  be  five  times 
faster.  However,  some  R&D  and  cost  overrun  problems  may  delay  its 
introduction.  Recently,  Cray’s  older  systems  have  begun  to  be  challenged  by 
comparatively  low-cost  mini-supercomputers  such  as  Convex  (of  Texas), 
Elxsi  (of  California),  Multiflow  Computer  Corp.  (of  Connecticut),  and 
Floating  Point  Systems  (of  Oregon).  These  machines  cost  approximately 
half  a  million  dollars  while  the  Cray’s  new  Y-MP  is  in  the  $20  million  range. 
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Many  predict  that  the  number  of  mini-supercomputers  will  grow  to  five 
times  between  1987  to  1992.  But  others  feel  that  the  market  is  already 
saturated  as  of  1988.  Figure  3  illustrates  the  computational  speed  of  major 
main  frame  computers  over  the  last  three  decades  while  Fig.  4  gives  a 
general  idea  of  the  correlation  between  the  computer’s  speed  and  the 
hardware  cost. 

According  to  the  designer  of  Cray-2,  Steve  Chen  (personal 
communication,  1988),  the  next  generation  of  vector  machines  which  is  being 
developed  by  his  IBM-supported  Super  Computer  Corporation  will  have  a 
processing  speed  of  nearly  one  hundred  times  the  Y-MP.  The  targeted  date 
for  introduction  will  be  approximately  five  years  from  now. 


VECTOR  PROGRAMMING  AND  THE  VECTORIZATION  OF  CFD 
CODES 

Even  though  some  compilers  can  perform  automatic  vectorization,  the 
final  processing  speed  of  most  CFD  codes  depends  on  a  programmer’s 
knowledge  of  the  machine’s  architecture,  so  that  optimum  concurrency  can 
be  obtained.  Numerical  schemes  used  in  CFD  modeling  also  influence  its 
processing  speed.  In  many  cases,  explicit  time  integration  schemes  vectorize 
naturally  because  they  are  an  algebraic  operation.  Implicit  numerical 
integration  schemes  require  additional  preparation  since  matrix  calculations 
are  involved. 

Most  vector  processors  contain  vector  pipelines.  Cyber-205  contains 
either  one,  two,  or  four  vector  pipelines,  each  of  which  is  fed  from  a  vector 
stream  unit.  The  top  speed  of  the  four-pipe  unit  is  400  mega-FLOPS  after 
the  initial  start-up  period. 

In  CFD  modeling,  the  objective  is  to  write  and  compile  the  numerical 
code  to  suit  the  machine’s  architecture.  Therefore,  the  modeler  should  have 
detailed  knowledge  about  the  hardware  system  as  well  as  aerodynamics  and 
aircraft  design  criteria.  Most  CFD  simulations  are  carried  out  by 
supercomputers  equipped  with  either  an  array  or  vector  processor. 

Unlike  a  vector  processor,  which  is  either  an  integral  part  of  the  CPU  or 
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resides  on  the  same  data  bus,  array  processors  coimect  to  the  I/O  bus  of  a 
computer.  The  speed  of  an  I/O  bus  can  be  an  order  of  magnitude  slower 
than  the  CPU  bus.  So  the  simulation  speed  is  much  slower  on  an  array 
processor  unless  the  array  is  extremely  large,  which  is  usually  the  case  for 
PNS  or  FNS  codes.  Array  processors  must  also  pay  an  additional  speed 
penalty  for  the  initialization  step  as  well  as  having  to  interact  with  the  host 
CPU. 

There  are  basically  two  ways  to  design  a  vector  processor: 
memory-to-memory  architecture  or  the  load/store  architecture.  The 
Control  Data  family  of  supercomputers  uses  the  memory-to-memory  method 
which  loads  data  directly  into  arithmetic  units  on  the  vector  board  from 
main  memory  without  the  need  for  vector  registers  or  the  so-called  cache. 
The  advantage  of  this  architecture  lies  in  its  suitability  for  very  large 
datasets. 

One  way  of  implementing  the  load/store  architecture  is  called  the 
"synchronous  method,"  which  is  adapted  by  Cray  and  IBM  3090.  In  this 
method,  a  cache  is  used.  During  CFD  simulation,  calculations  are  carried 
out  following  the  main,  scalar  CPU  until  a  vector  call  is  encountered.  While 
the  vector  is  being  processed,  the  CPU  is  in  a  wait  state.  The  CPU’s  wasted 
idle  time  is  directly  proportional  to  the  size  of  the  vector.  The  newest 
machines  from  Convex  and  DEC  (C240  and  VAX  9000,  around  the  spring  of 
1990)  have  adapted  the  asychronous  load/Store  method.  This  type  of 
architecture  allows  the  main  CPU  to  perform  in  paraUel  with  the  vector 
processor  if  the  operations  in  the  CPU  do  not  depend  on  the  results  from  the 
vector  processor. 

In  terms  of  programming  software,  FORTRAN  is  still  the  most  popular 
scientific  programming  language.  There  are  special  utility  programs  used  by 
many  super-  and  mini-supercomputer  manufacturers,  such  as  ETA,  and 
CYBER  (CDC),  to  vectorize  FORTRAN  programs  to  take  advantage  of  the 
vector  processing  hardware.  One  vectorization  utility  is  called  VAST 
(vector  and  array  syntax  translator)  and  is  produced  by  Pacific  Sierra 
Research.  To  give  a  simple  example,  if  a  constant  is  to  be  added  to  every 
element  of  an  array  using  scalar  FORTRAN,  th.  ‘llovring  DO  loop  is 
needed. 
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DO  10  1=1,2000 
10  A(I)  =  A(I)  -f  X 

When  this  operation  is  performed  on  an  ETAIO  which  contains  hardware 
specially  designed  for  vector  operation,  the  ETA  VAST-2  can  transform  the 
above  scalar  operation  into; 

A(  1,2000)  =  A(  1,2000)  -b  X 

The  FORTRAN  compiler  of  ETAIO  would  recognize  this  vector 
operation  and  will  generate  the  appropriate  machine  instructions  to  load  the 
vector  processor  and  add  the  constant  to  all  elements  of  the  array 
simultaneously.  The  vectorized  code  would  take  about  only  l/2000th  the 
time  to  complete  the  operation. 

Since  vectorization  can  be  accomplished  in  several  ways  for  a  complex 
program,  the  preprocessor  must  determine  the  most  efficient  one.  Even 
though  the  vectorization  preprocessor  can  be  a  powerful  tool,  it  sometimes 
can  give  unwanted  vectorizations  resulting  in  harmful  side  effects. 
Programmers  sometimes  have  to  give  directives  in  the  source  file  that  help 
syntex  translators  such  as  VAST  make  the  vectorization. 


CFD  MODELING  ON  ARRAY  COMPUTERS  AND  CONNECTION 
MACHINES 

An  array  computer  consists  of  many  microprocessors  that  perform 
identical  computations.  The  concept  is  originated  from  von  Neumann. 
ILLIAC  IV,  which  consists  of  64  (8x8)  processing  elements,  is  the  earliest 
example  of  a  general-purpose  array  computer.  Even  though  each  processor 
in  the  ILLIAC  IV  is  quite  powerful,  the  size  of  the  array  is  too  small  to 
perform  practical  calculations  in  CFD  modeling.  In  1976,  RAND  proposed 
to  build  a  Navier  Stokes  array  computer  consisting  of  10,000  identical 
microprocessors  arranged  in  an  100x100  array  (Gritton  et  al.,  1977).  The 
concept  was  to  perform  CFD  simulations  on  an  array  computer  that  is 
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designed  only  to  solve  the  Navier  Stokes  equation  with  great  efficiency. 

They  proposed  to  solve  the  three-dimensional  Navier  Stokes  equations  in 
primitive  variables  and  the  Poisson  equation  for  pressure  via  Fourier 
transformation  in  the  span-wise  direction.  With  100x100  array  of 
microprocessors  and  128  Fourier  components,  the  processing  speed  was 
expected  to  be  several  orders  faster  than  a  sequential  machine.  The  main 
reason  for  this  difference  is  that  during  a  simulation,  most  of  the  hardware  in 
the  sequential  processor  sits  idle.  In  a  parallel  machine,  calculations  are 
carried  out  at  every  processor  simutaneously. 

Several  general-purpose,  massively  parallel  processors  such  as  the  one 

built  at  the  Thinking  Machine  Corp.  (Hillis,  1985)  are  presently  being 

developed  and  used  for  CFD  modeling.  For  example,  the  Naval  Research 

Laboratory  has  two  machines  presently  dedicated  to  CFD  simulations.  The 

smaller  one  (nicknamed  Bambi)  consists  of  8000  floating-point  processing 

chips  in  a  4K-I-4K  congifuration.  The  larger  one  (nicknamed  Godzilla)  has 

16,000  floating-point  processing  chips  in  a  8K-J-8K  configuration.  Each 

processor  has  2K  words  of  memory.  Navier  Stokes  solutions  coded  on  these 

* 

massively  parallel  connection  machines  are  usually  written  in  C  language, 
which  is  quite  similar  to  FORTRAN.  Usually,  C  language  is  used  together 
with  PARIS  or  CMIS  commands  to  increase  the  processing  speed.  Large 
connection  machines  can  sometime  outperform  supercomputers.  For 
example,  a  256xl28xl28-grid,  three-dimensional  fluid  simulation  of  complex 
shock  interaction,  when  carried  out  using  the  16K  Connection  Machine  at 
the  Naval  Research  Laboratory,  takes  20  seconds  per  time  step.  The 
identical  problem  would  take  80  seconds  per  time  step  on  a  Cray  X-MP 
(Boris  et  al.,  1989).  This  implies  that  the  speed  of  the  16K  connection 
machine  is  equivalent  to  the  speed  of  Cray’s  Y-MP. 

Recent  computations  of  fully  compressible  two-dimensional  Navier 
Stokes  equations  (1024  x  2048  grid)  on  the  large  parallel  connection  machine 
CM2  at  the  University  of  Colorado  achieved  a  speed  of  1,500  Mega-FLOPS 
(1.5  G-FLOPS,  Biringen  et  al.,  1989).  The  benchmark  test  was  carried  out 
using  the  MacCormack  second-order  explicit  scheme  for  the  simulation  of  a 
square  wave  propagating  outward  through  a  square  domain  with  radiation 
boundary  conditions  for  24,000  integration  time  steps.  Coded  in  C  and 
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PARIS,  the  execution  speed  in  the  parallel  connection  machines  slows  down 
substantially  when  implicit  schemes  are  used  because  of  the  involvement  of 
matrix  operations. 

Presently,  on  connection  machines,  the  finite-difference  scheme  is  more 
efficient  than  the  spectral  method  (Sec.  3)  for  lower  accuracy  modeling  such 
as  second-  or  fourth-order  accuracy.  When  the  required  accuracy  increases, 
the  pseudospectral  method  is  better,  particularly  for  simpler  geometric 
shapes.  Recently  many  CFD  simulations  of  the  hydrodynamics  of 
symmetrical  and  spherical  systems  have  been  studied  using  parallel  spectral 
methods  on  connection  machines  (Pelz,  1989). 
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7.  ASPECTS  OF  MODELING  NEEDS  AND  UNCERTAINTIES  IN 
HYPERSONIC  SIMULATION 

REQUIREMENTS  ON  THE  COMPUTATIONAL  SPEED  AND  MEMORY 

CFD  models  are  based  on  the  approximated  (numerical)  solutions  of  the 
governing  equations  over  the  finite-grid  network.  Consequently,  the  finer 
the  grid,  the  more  accurate  the  numerical  solutions  are.  There  are  tradeoff 
between  speed  and  computer  memory  requirements.  To  obtain  the  best 
estimate  of  the  flow  field  around  an  aircraft,  we  often  trade  computational 
speed  for  higher  spatial  resolution  with  the  largest  machine  available.  To  do 
this,  we  sometimes  have  to  use  a  less  efficient  but  more  straightforward 
numerical  method  such  as  the  explicit  scheme.  This  type  of  scheme  requires 
fewer  variables  residing  within  the  machine’s  main  memory  at  the  same 
time.  Variables  that  are  not  involved  reside  in  the  machine’s  external 
memory  units  waiting  to  be  called  upon.  More  "time-efficient"  and 
sophisticated  integration  methods  such  as  the  implicit  method  would  require 
more  variables  in  the  main  memory  simultaneously.  At  the  present  time, 
the  largest  computer’s  main  memory  and  speed  can  allow  simulations  with  a 
half-million  to  a  million  grid  points  if  the  simplest  expUcit  scheme  is  used. 

In  this  type  of  scheme,  only  five  neighboring  grid  points  are  contained  in  the 
CPU  at  a  given  time  during  the  time  and  spatial  marching  (integration) 
process.  A  paging  process  is  involved  in  which  cyclic  data  blocks  move  in 
and  out  of  the  CPU  and  the  external  high  speed  storage  device. 

With  the  maximum  allowable  resolution  using  the  method  described 
above,  the  accuracy  of  the  CFD  approach  at  lower  Mach  range  is  around  5  to 
7  percent.  For  example,  in  a  simulation  involving  an  X24C-10D 
experimental  plane,  the  estimated  values  of  the  lift  and  drag  coefficients, 
using  the  CFD  method  as  compared  vnth  experiments  at  M=5.95,  are  4.71 
percent  and  6.71  percent,  respectively  (Shang  and  Scherr,  1985).  This 
simulation  was  made  vrith  the  parabolized  Navier  Stokes  equation  coded 
with  an  explicit  scheme  (MacCormack,  1969)  coupled  with  a  zeroth-order 
algebraic  turbulence  closure  technique  (Baldwin  and  Lomax,  1978). 
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When  the  computed  values  are  compared  with  the  measurements,  the 
larger  deviations  are  located  near  the  sharper  geometric  transitions.  This 
indicates  that  higher  spatial  resolution  will  improve  the  accuracy.  A 
higher-order  turbulence  model  will  also  improve  results  around  curved 
surfaces.  All  of  this  indicates  a  need  for  bigger  and  faster  computers. 

Within  the  foreseeable  future,  the  computational  speed  of  the  computer 
is  likely  to  follow  the  previous  trend,  namely,  to  increase  3  to  5  times  in  3  to 
4  years.  The  present  computational  speed  limit  is  around  1  billion  floating 
point  operations  per  second  (FLOPS).  In  projecting  future  speeds  of 
supercomputers,  the  improvements  can  be  made  in  several  ways.  In  terms  of 
hardware  development,  the  speed  improvement  due  to  the  increase  in  chip 
density  is  expected  to  increase  by  tenfold  every  seven  years  (Fig.  3). 
Supercomputers  can  also  increase  their  processing  speed  through  parallelism 
in  architectural  design,  and  the  I/O  speed  can  be  increased  by  the  use  of 
super-conductive  material,  just  to  name  a  few  possibibties. 

For  the  NASP  and  the  N ASP-derived  system  applications  around  the 
year  2000,  the  expected  computational  speed  is  in  the  neighborhood  of  10-20 
billion  FLOPS.  A  system  with  this  capability  could  carry  out  CFD 
simulations  for  an  external  aircraft  configuration  consisting  of  10  million  grid 
points  assuming  the  use  of  numerical  techniques  similar  to  those  described 
above. 


REQUIREMENTS  ON  THE  SPECIFICATION  OF  BOUNDARY 
CONDITIONS 

The  use  of  CFD  is  based  on  the  solution  of  boundary  value  problems.  To 
obtain  a  solution,  boundary  conditions  have  to  be  specified.  Any  inaccuracy 
in  specifying  the  proper  boundary  condition  propagates  into  the 
computational  field.  But  if  one  knows  the  exact  conditions  at  the  boundary, 
he  would  already  have  the  solution  to  the  problem.  This  is  a  well-known 
difficulty  in  the  field  of  numerical  modeling.  One  solution  would  be  to  set 
the  boundaries  of  the  model  as  far  as  possible  from  the  aircraft  so  that  its 
influence  will  be  very  small  near  the  far-field  boundaries.  This  approach 
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would  need  a  large  computational  capacity.  The  next  solution  is  to  design  a 
sophisticated  numerical  scheme  at  the  model’s  boundaries  in  such  a  way  that 
the  dynamic  effects  (e.g.,  shock)  originated  at  the  aircraft  would  leave  the 
boundaries  without  reflecting  back  into  the  model  domain.  This  type  of 
boundary  treatment  is  called  the  "radiation  boundary  condition."  It  means 
that  the  energy  within  the  model  can  radiate  out  of  the  model  through  its 
boundaries. 

To  make  the  boundary  condition  completely  radiative  is  easier  said  than 
done.  One  major  difficulty  is  to  design  a  numerical  scheme  at  the  boundary 
which  can  "sense"  all  outgoing  (pressure)  waves  of  various  frequencies  and 
let  them  pass  through.  Furthermore,  it  is  even  more  difficult  to  design  such 
a  scheme  if  "zonal  modeling"  techniques  are  used.  In  such  a  method,  several 
models  are  patched  together  to  form  an  entire  aircraft.  Mutual  boundaries 
are  located  throughout  the  system  and  all  of  them  need  to  be  treated.  Since 
the  boundaries  interconnect,  special  numerical  schemes  are  therefore  needed 
to  overcome  the  redundancy  in  their  specification.  This  is  also  true  in  the 
specification  of  the  radiation  boundary  conditions  for  an  adaptive  (moving) 
grid  network.  To  avoid  this  type  of  problem,  it  is  much  easier  to  simulate 
the  entire  aircraft  in  a  single  model,  if  the  computational  resources  are 
available. 


REQUIREMENTS  ON  THE  VERIFICATION  DATA  AT  THE  UPPER 
HYPERSONIC  RANGE 

During  different  stages  of  model  development,  experimental  data  are 
needed  to  verify  the  theoretical  basis  and  the  computational  code  of  the 
numerical  model.  As  the  Reynolds  number  and  air  speeds  increase,  the 
accuracy  of  the  experimental  measurement  decreases.  The  degree  of 
difficulty  of  acquiring  the  experimental  data  also  increases  as  the  sj>eed 
increases.  Shock  tubes  and  gun-launched  models  can  reach  high  Mach 
range,  but  they  all  have  limitations  for  NASP  applications.  Because  of  the 
lack  of  emphasis  on  the  hypersonic  experiments  during  the  1970s,  hypersonic 
wind  tunnels  and  other  ground  test  facilities  designed  to  perform  continuous 
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free  flight  tests  are  either  no  longer  in  existence  or  out  of  commission.  At 
present,  very  few  ground  testing  facilities  or  wind  tunnels  can  simulate 
continuous  free  flight  conditions  at  Mach  number  higher  than  8.  At  speeds 
above  Mach  10,  it  becomes  increasingly  more  difficult  to  obtain  high  enough 
temperature  to  avoid  gas  condensation  with  corresponding  high  pressure  to 
simulate  flight  Reynolds  number  range.  Recently,  tetrafluoromethane  (CF^) 

and  nitrogen  (N2)  have  been  used  (Midden  and  Miller,  1978)  for  simulating 

the  thin  shock  layer,  chemically  reactive,  and  high  Reynolds  number  effects 
associated  with  the  hypersonic  flows. 


REQUIREMENTS  ON  THE  UNIVERSALITY  OF  A  MODEL'S 
PREDICTION  CONSTANTS 

Without  suitable  verification  data  at  high  speed,  the  predictabiUty  or 
accuracy  of  the  CFD  codes  will  depend  upon  the  universality  of  their 
turbulence  model.  In  other  words,  the  values  of  the  closure  constants  have 
to  stay  the  same  for  these  models  to  predict  beyond  the  verified  range  of 
applicability.  This  seems  to  be  one  of  the  crucial  issues  at  the  present  time. 
New  CFD  codes  are  undoubtedly  versatile  and  powerful.  At  low  Mach 
range,  experimental  results  have  been  used  extensively  to  verify  and  adjust 
turbulence  closure  computations.  Since  turbulence  modeling  is  not 
system-specific  in  terms  of  geometric  shape,  it  has  been  a  favored  method  for 
model  verification.  During  the  adjustment  process,  the  so-called  turbulence 
closure  "constants"  have  been  optimally  fitted  for  various  proposed  models. 
Some  of  the  results,  presented  in  Table  5,  show  some  degree  of  uncertainty. 
At  this  moment,  we  cannot  predict  whether  the  degree  of  uncertainty  will 
increase  in  the  higher  Mach  range.  To  obtain  verification  data  at  high  speed 
(M  8-25),  measurements  from  shuttle  reentry  flight  or  drones  launched  from 
high-speed  planes  may  provide  some  technical  data  for  benchmarking 
computational  codes.  In  turbulence  modeling,  areas  of  low  predictive 
reliability  are  listed  in  Tables  6  and  7. 

In  the  field  of  fluid  mechnics,  however,  there  are  very  few  truly 
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"universal"  constants.  Constants  often  considered  as  universal,  such  as  the 
Kolmogorov  constant  are  later  found  to  vary  under  different  flow  conditions. 
The  Kolmogorov  constants  for  the  distribution  of  turbulent  spectral  energy 
within  turbulent  flow  was  recently  found  to  vary  in  two-dimensional 
turbulence  (Kraichnan,  1987;  Qian,  1986).  This  finding  is  quite  important 
since  flows  under  the  thin-layer  approximation  or  flows  describable  by 
parabolized  CFDs  are  two-dimensional  turbulences.  To  simulate  the  proper 
laminar/turbulent  transition  process  at  high  hypersonic  flight  speeds,  models 
of  compressible  turbulence  should  be  used.  For  example,  the  entrainment 
rate  of  compressible  turbulence  is  lower  than  the  corresponding 
incompressible  turbulence.  The  effects  of  shock  structure  on  the  transport  of 
kinetic  energy  in  and  out  of  a  turbulent  eddy  is  still  not  completely 
understood  (Breidenthal,  1990).  Furthermore,  the  compressibility  of  air  at 
hypersonic  speed  may  change  the  universally  accepted  concept  of 
Kolmogorov’s  law  of  the  cascade  of  spectral  energy  in  the  frequency  domain 
(Dimotakis,  1989).  It  is  still  too  early  to  know  the  impact  of  these  new 
findings  on  CFD  modeling. 


REQUIREMENTS  ON  THE  COOPERATION  BETWEEN 
MODELERS  AND  DESIGNERS 

To  develop  and  work  with  CFD  methods  and  codes,  one  needs  an 
extensive  background  in  both  mathematics  and  computer  sciences  as  well  as 
fluid  dynamics.  Sensitivity  to  experimental  findings  is  also  essential.  CFD 
modeling  gradually  evolves  into  a  highly  specialized  field.  To  maintain  a 
CFD  simulation  system  often  becomes  a  full-time  task  for  a  small  research 
team.  After  completion  of  the  Apollo  moon-landing  program  and  design  of 
the  shuttle  reentry  vehicle  in  the  early  1970s,  there  has  been  no  major  design 
activity  in  hypersonics  during  the  last  one  and  half  decades.  As  a 
consequence,  most  CFD  modelers  have  limited  design  experience.  Since 
computers  are  not  yet  fast  and  big  enough  to  perform  a  complete  automated 
design  under  given  constraints,  CFD  simulations  still  play  the  role  of  the 
designer’s  aids,  which  are  extremely  complicated  to  use.  To  close  this  gap. 
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experienred  designers  have  to  work  with  the  CFD  simulation  team.  This  is 
necessary  because  an  experienced  designer  not  only  can  give  valuable 
guidance  in  the  design  phase  but  can  also  check  to  see  if  the  computer 
output  makes  sense.  As  mentioned  earlier,  although  they  are  powerful  tools, 
CFD  results  are  still  only  simulations. 


COMPATIBILITY  REQUIREMENTS  BETWEEN 
HARDWARE  AND  SOFTWARE 

The  computational  efficiency  of  a  numerical  scheme  is  a  relative  matter. 
It  often  depends  on  the  type  of  hardware  the  code  is  running  on.  An 
understanding  of  the  machine  architecture  is  very  important  for  code 
development.  For  example,  during  the  numerical  integration  process,  the 
control  of  data  flow  between  the  parallel  functional  unit  and  the 
hierarchically  organized  memory  is  essential.  Because  the  computer’s  basic 
clock  speed,  which  determines  the  rate  of  calculation,  varies  from  machine  to 
machine,  the  cycle  time  required  for  various  computer  functions  will  also 
vary  thereby  resulting  in  major  differences  in  overall  computational  speed. 
Machine-specific  characteristics  such  as  the  "memory  cycle  time,"  the 
"access  time,"  and  the  "transfer  rate"  all  require  a  different  number  of  clock 
periods  to  complete.  A  CFD  code  should  be  designed  to  allow  maximum 
concurrent  operations  to  achieve  the  overall  efficiency. 

According  to  Cray  Research  (Shang  et  al.,  1980),  the  efficiency  of  a 
FORTRAN  code  can  be  improved  by: 

•  Avoiding  a  long  and  complicated  "do  loop" 

'  Replacing  the  most  frequently  addressed  temporary  scalar  variable 
with  a  vector  temporary  variable 

•  Developing  the  code  to  allow  maximum  use  of  concurrent 
(chaining)  operations.  The  chaining  of  two  or  more  vector 
operations  allows  a  vector  machine  to  yield  more  than  one  result 
per  clock  period. 
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Along  the  same  line  as  discussed  above  but  more  fundamental  in  nature, 
it  is  important  to  use  better  programming  languages  and  operating  systems. 
At  present,  the  majority  of  aerospace  engineers  write  CFD  models  in 
FORTRAN,  which  is  considered  to  be  primitive  as  compared  to  the  more 
advanced  languages  such  as  Pascal  or  ADA.  More  efficient  transfer  between 
the  database  and  the  arithmetic  unit  can  also  be  achieved  by  improving  the 
operating  system. 
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8.  CONCLUSIONS 


During  this  review  the  following  conclusions  were  reached: 

1.  Computational  fluid  dynamic  (CFD)  models,  if  they  are  properly 
designed,  verified,  and  applied,  are  versatile  and  powerful  tools.  With  the 
level  of  resolution  that  can  be  handled  by  the  present  state-of-the-art 
computers,  the  level  of  verifiable  simulation  accuracy  (in  terms  only  of 
pressure  distribution)  is  approximately  within  5  to  6  percent  at  the 
intermediate  speed  range  (Mach  0  -  8).  But,  it  is  difiBcult  to  even  estimate 
the  inaccuracies  in  such  crucial  parameters  as  skin  frictions,  heating  values, 
and  turbulent  mixing.  When  the  computed  values  are  compared  with  the 
available  data,  the  larger  deviations  are  located  near  the  sharper  geometric 
transitions.  This  indicates  that  higher  spatial  resolution  will  improve  the 
accuracy.  A  higher-order  turbulence  model  should  also  improve  results 
around  curved  surfaces.  These  improvements  need  bigger  and  faster 
computers.  However,  it  may  be  years  before  we  can  narrow  uncertainties  in 
both  the  laminar-turbulent  transition  location  and  the  characteristics  of 
flows  as  they  enter  the  inlet  region.  We  may  also  require  time  to  better 
understand  the  character  of  the  mixing  region  with  the  SCRAMJET  as  the 
Mach  numbers  approach  those  of  interest  to  NASP. 

2.  For  better  predictive  accuracy  and  for  the  generality  in  the  CFD 
models’  predictive  applicability,  turbulence  closure  techniques  need  to  be 
further  improved.  Experimental  data  at  the  full  range  of  NASP  speeds  are 
needed  to  evaluate  the  closure  constants.  A  wider  range  of  applicability  can 
be  achieved  if  the  constants  are  more  universal.  Universality  implies  higher 
predictive  reliability  into  higher  speed  range.  Closure  models  will  nearly 
always  be  needed  in  CFD  within  the  foreseeable  future  (e.g.,  until  2000) 
unless  a  major  breakthrough  in  computing  power  takes  place.  The  speed  and 
memory  of  computers  would  have  to  inprove  several  orders  of  magnitude 
over  present  state-of-the-art  machines  considering  realistic  flight  Reynolds 
number  range.  With  improvement  of  this  magnitude,  the  troublesome 


-101- 


"turbulence  closure"  step  can  then  be  completely  eliminated.  Given  that  the 
past  rate  of  progress  in  computing  power  continues  (i.e.,  10  times  in  7  years) 
and  the  NASP  schedule  remains  unchanged,  it  will  be  necessary  to  continue 
using  turbulence  modeling  for  the  entire  NASP  development  program,  as 
well  as  possibly  for  the  NDV  development  program  as  well.  To  approach 
universality  for  constants,  many  more  data  must  be  collected  at  all  NASP 
speed  ranges.  If  this  is  not  done,  a  similar  level  of  predictive  uncertainty  as 
the  present  will  continue  to  exist. 

3.  Because  of  the  relatively  inactive  period  in  hypersonic  research  for  the 
past  15  years,  the  design-supporting  ground-testing  facilities  that  were  in 
existence  have  not  been  properly  maintained  or  are  out  of  commission.  At 
present,  there  is  an  urgent  need  for  ground-testing  facilities  with  improved 
measurement  techniques  to  provide  verification  data  for  CFD  models. 
Complementary  techniques  to  ground  tests  are  needed  to  verify  CFD 
simulations  at  the  high  Mach  ranges  where  experimental  data  are  difficult  to 
obtain.  Recent  advances  in  flow  measurements  by  means  of  laser-Doppler 
technology  will  undoubtedly  improve  the  quality  of  verification  data  in  the 
future.  The  shock  tube  "RHYFL"  (Rockdyne  Hypersonic  Facility),  which 
will  be  completed  in  the  near  future,  has  a  velocity  limit  of  Mach  16.  The 
piston  is  driven  by  nitrogen  with  the  testing  duration  of  1/1000  second.  As  a 
consequence,  dynamic  processes  induced  by  transient  dynamics  may  pose 
some  problem.  But  extensive  use  of  RHYFL  will  be  an  important  part  of 
reducing  the  uncertainty  in  the  CFD  simulation. 

4.  The  numerical  solution  schemes  used  in  the  CFD  models  have 
improved  in  accuracy  and  efficiency  during  the  last  two  decades.  However, 
no  clear  distinction  has  developed  among  all  the  available  methods  as  to 
which  is  the  best  overall  method.  Tradeoffs  often  have  to  be  made.  There  are 
times,  when  less  efficient  but  more  straightforward  integration  schemes  have 
been  selected  so  that  fewer  spatial  points  reside  in  the  CPU  simultaneously. 
By  doing  this,  a  great  number  of  grid  points  can  be  simulated  with  a  single 
model,  thus  simplifying  the  boundary  condition  handling  problem,  even 
though  this  approach  increases  the  simulation  time  (I/O  process). 

5.  One  major  advancement  in  the  field  of  CFD  is  the  application  of 
coordinate  transformation  schemes  often  employed  in  conjunction  with  the 
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finite-difference  methods.  When  an  implicit  scheme  is  used  together  with 
transformed  coordinates,  the  drop  in  accuracy  of  numerical  approximation 
should  be  evaluated  especially  near  the  crucial  area  such  as  the  leading  edge 
where  the  Courant  number  is  very  high  as  a  result  of  small  grid  size. 

6.  During  the  coming  years,  for  NASP  applications,  emphasis  on  CFD 
research  should  be  on  the  following  crucial  areas; 

Validation  of  higher  order  turbulence  models  to  reduce  the 
level  of  uncertainty  in  the  areas  of  low  predictive  reliability  such  as 
rapid  compression/expansion,  strong  swirl,  kinematically  influenced 
chemical  reaction,  dynamic  instability. 

Improvement  in  modeling  compressible  turbulence. 

Improvement  in  predicting  the  location  and  the  length  of  the 
boundary  layer  transition  zone  on  the  aerospace  plane  and  the  details 
of  the  transition  process. 

Improvement  in  conserving  mass,  momentum,  and  energy  in  the 
numerical  scheme  solving  the  Navier  Stokes  equations  involving  grid 
transformation  particularly  curvilinear  transformation. 

Treatment  of  boundary  conditions  particularly  in  a  "nested" 
modeling  system  where  different  grids  are  used. 

•  More  efficient  cooperation  between  government  and  industry  in  CFD 
code  development  and  validation. 

Explicitly  evaluating  and  publishing  the  uncertainty  in  CFD 
simulation  results  as  a  function  of  the  vehicle  speed  and  position 
along  the  vehicle. 
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